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ABSTRACT 

We exploit the deep and extended far-IR data sets (at 70, 100 and 160 \im) of the Herschel 
Guaranteed Time Observation (GTO) PACS Evolutionary Probe (PEP) Survey, in combination 
with the Herschel Multi-tiered Extragalactic Survey data at 250, 350 and 500 pm, to derive the 
evolution of the rest-frame 35-, 60-, 90- and total infrared (IR) luminosity functions (LFs) up 
to z ~ 4. We detect very strong luminosity evolution for the total IR LF (Lir a (1 + z) 3 ' 55 dt 010 
up to z ~ 2, and a (1 + z) L62±0 ' 51 at 2 < z < 4) combined with a density evolution (a 
(1 + o.57 ± 0.22 U p to z ~ 1 and oc (1 + z )-3.92±o.34 a t l < z < 4). i n agreement with 

previous findings, the IR luminosity density (pir) increases steeply to z ~ 1, then flattens 
between z ~ 1 and z ~ 3 to decrease at z > 3. Galaxies with different spectral energy 
distributions, masses and specific star formation rates (SFRs) evolve in very different ways 
and this large and deep statistical sample is the first one allowing us to separately study the 
different evolutionary behaviours of the individual IR populations contributing to pir. Galaxies 
occupying the well-established SFR-stellar mass main sequence (MS) are found to dominate 
both the total IR LF and pir at all redshifts, with the contribution from off-MS sources (>0.6 
dex above MS) being nearly constant (~20 per cent of the total pir) and showing no significant 
signs of increase with increasing z over the whole 0.8 < z < 2.2 range. Sources with mass in 
the range 10 < log(M/Mo) <11 are found to dominate the total IR LF, with more massive 
galaxies prevailing at the bright end of the high-z (>2) LF. A two-fold evolutionary scheme 
for IR galaxies is envisaged: on the one hand, a starburst-dominated phase in which the 
Super Massive Black Holes (SMBH) grows and is obscured by dust (possibly triggered by a 
major merging event), is followed by an AGN-dominated phase, then evolving towards a local 
elliptical. On the other hand, moderately star-forming galaxies containing a low-luminosity 
AGN have various properties suggesting they are good candidates for systems in a transition 
phase preceding the formation of steady spiral galaxies. 

Key words: galaxies: active -galaxies: evolution -galaxies: luminosity function, mass func- 
tion-galaxies: starburst- cosmology: observations -infrared: galaxies. 
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1 INTRODUCTION 

Understanding the origin and growth of the galaxies we observe 
today is one of the main problems of current cosmology. The lu- 
minosity function (LF) provides one of the fundamental tools to 
probe the distribution of galaxies over cosmological time, since it 
allows us to assess the statistical nature of galaxy formation and 
evolution. When computed at different redshifts, the LF constitutes 
the most direct method for exploring the evolution of a galaxy 
population, describing the relative number of sources of different 
luminosities counted in representative volumes of the Universe. The 
LF computed for different samples of galaxies can provide a crucial 
comparison between the distribution of different galaxy types, i.e. 
galaxies at different redshifts, in different environments or selected 
at different wavelengths. 

It has now become clear that we cannot understand galaxy evo- 
lution without accounting for the energy absorbed by dust and re- 
emitted at longer wavelengths (i.e. Genzel & Cesarsky 2000), in 
the infrared (IR) or sub-millimetre (sub-mm). Dust is responsible 
for obscuring the ultraviolet (UV) and optical light from galaxies: 
since star formation (SF) occurs within dusty molecular clouds, far- 
IR and sub-mm data, where the absorbed radiation is re-emitted, 
are essential for providing a complete picture of the history of SF 
through cosmic time, which is one of the fundamental instruments 
we have to reconstruct how galaxies have evolved since their forma- 
tion epoch. For these reasons, extragalactic surveys in the rest-frame 
IR represent a key tool for understanding galaxy formation and 
evolution. 

Surveys of dust emission performed with the former satellites 
exploring the Universe in the mid- and far-IR domain, i.e. the In- 
frared Astronomical Satellite (IRAS; Neugebauer et al. 1984) and 
the Infrared Space Observatory (Kessler et al. 1996), allowed the 
first studies of the IR-galaxy LF at z < 0.3 (Saunders et al. 1990) 
and z < 1 (Pozzi et al. 2004), respectively. With Spitzer 24- pan 
data, it was possible to study the evolution of the mid-IR LF up to 
z ~ 2 (e.g. Le Floc’h et al. 2005; Caputi et al. 2007; Rodighiero 
et al. 2010a), while, even with the deepest Spitzer Space Telescope 
(Werner et al. 2004) 70- pm data, only z ~ 1-1 .2 could be reached in 
the far-IR (Magnelli et al. 2009; Patel et al. 2013) - though Magnelli 
et al. (2011) reached z ~ 2 through stacking. Since the rest-frame IR 
spectral energy distributions (SEDs) of star-forming galaxies and 
AGN peak at 60-200 pm, to measure their bolometric luminos- 
ity and evolution with z we need to observe in the far-IR/sub-mm 
regime. However, the detection of large numbers of high-z sources 
at the peak of their IR SED was not achievable before the Herschel 
Space Observatory (Pilbratt et al. 2010), due to source confusion 
and/or low detector sensitivity, and our knowledge of the far-IR LF 
in the distant Universe is still affected by substantial uncertainties. 
Ground-based and balloon-borne observations in the mm/sub-mm 
range, probing the evolution of the most distant (z > 2) and lumi- 
nous dusty galaxies, have so far been limited to the identification of 
sources at the very bright end of the LF (e.g. Chapman et al. 2005). 
All of these works detected strong evolution in both luminosity 
and/or density, indicating that IR galaxies were more luminous 
and/or more numerous in the past. Strong observational evidence 
of high rates of evolution for IR galaxies has been obtained also 
through the detection of a large amount of energy contained in the 
Cosmic Infrared Background (CIRB; Hauser & Dwek 2001), and 
the source counts from several deep cosmological surveys (from 15 
to 850 pm) largely exceeding the no-evolution expectations (e.g. 
Smail, Ivison & Blain 1997; Elbaz et al. 1999; Gruppioni et al. 2002; 
Papovich et al. 2004; Bethermin et al. 2010; Marsden et al. 2011). 


Both the CIRB and the source counts require a strong increase in 
the IR energy density between the present time and z ~ 1-2. At 
higher redshifts the total emissivity of IR galaxies is poorly con- 
strained, due to the scarcity of Spitzer galaxies at z > 2, the large 
spectral extrapolations to derive the total IR luminosity from the 
mid-IR (see e.g. Elbaz et al. 2010, Nordon et al. 2010, 2012 for 
descriptions of the failure, at least at z > 1.5, of previous total IR 
luminosity extrapolations from the mid-IR, although we must note 
that this failure mainly affects luminosity-dependent methods like, 
e.g., that of Chary & Elbaz 2001) and the incomplete information 
on the z-distribution of sub-mm sources (Chapman et al. 2005). 

Herschel , with its 3.5-m mirror, is the first telescope which allows 
us to detect the far-IR population to high redshifts (z ~ 4-5) and to 
derive its rate of evolution through a detailed LF analysis. The new 
extragalactic surveys provided by Herschel in the far-IR/sub-mm 
domain, like the wide and shallow Herschel - Astrophysical Tera- 
herz Large Area Survey (ATLAS) (Eales et al. 2010; Dunne et al. 
201 1), the complementary Herschel Multi-tiered Extragalactic Sur- 
vey (HerMES; Oliver et al. 2012) and PACS Evolutionary Probe 
(PEP; Lutz et al. 2011) covering the most popular cosmological 
fields, and the deep, pencil beam, Herschel-G reat Observatories 
Origin Deep Survey (GOODS) project (Elbaz et al. 2011), will be 
crucial to assess galaxy and AGN evolution in the IR at z > 2. They 
will give us the opportunity to study in detail the population of IR 
galaxies and their evolution with cosmic time since the Universe 
was about a billion years old. In particular, the Photodetector Ar- 
ray Camera and Spectrometer (PACS; Poglitsch et al. 2010), with 
its high sensitivity and resolution at 70-, 100- and 160-qm, is the 
best- suited instrument to detect faint IR sources by overcoming the 
source confusion and blending problems that affected the previous 
far-IR missions. 

This is the first of two papers aiming at deriving the far-IR and 
total IR LFs from the Herschel PACS + Spectral and Photometric 
Imaging Receiver (SPIRE; Griffin et al. 2010) data obtained within 
the PEP and HerMES extragalactic survey projects. In this paper, 
we derive the rest-frame 35-, 60- , 90- and total IR (8-1000 jam) LFs 
from a sample selected at PACS 70, 100 and 160 pan wavelengths in 
the GOODS (GOODS-S and GOODS-N), Extended Chandra Deep 
Field South (ECDFS) and Cosmic Evolution Survey (COSMOS) 
areas. We use the full 70-500 jam PACS+SPIRE data to determine 
Lir and SED properties of the PACS-selected sources. In a related 
paper, Vaccari et al. (in preparation) derive rest-frame 100-, 160- and 
250- jam and total IR LFs for a SPIRE- selected sample. In addition, 
a third work aimed at studying the total IR LF based on the 24- jam- 
selected sample, using all the PEP+HerMES data in the COSMOS 
field, is ongoing (Le Floc’h et al., in preparation). 

PEP is one of the major Herschel Guaranteed Time extragalactic 
key-projects, designed specifically to determine the cosmic evolu- 
tion of dusty SF and of the IR LF. It is structured as a ‘wedding 
cake’, based on four different layers covering different areas to dif- 
ferent depths at 100 and 160 jam (in the GOODS-S field also at 
70 jam), from the large and shallow COSMOS field to the deep, 
pencil beam GOODS-S field. PEP includes the most popular and 
widely studied extragalactic fields with extensive multiwavelength 
coverage available, in particular deep optical, near-IR and Spitzer 
imaging and spectroscopic and photometric redshifts: COSMOS, 
Lockman Hole, Extended Groth Streep, ECDFS, GOODS-N and 
GOODS-S (see Berta et al. 2010, 2011 and Lutz et al. 2011 for 
a detailed description of the fields and observations). Coordinated 
observations of the PEP fields at 250, 350 and 500 jam with SPIRE 
have been obtained by the HerMES (Oliver et al. 2012). HerMES, 
analogously to PEP but extending to a much wider area, is a legacy 
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programme designed to map a set of nested fields (~380 deg 2 in 
total) of different sizes and depths, using SPIRE (at 250, 350 and 
500 jam), and PACS (at 100 and 160 jam, shallower than PEP), with 
the widest component of 270 deg 2 with SPIRE alone. In the fields 
covered by PEP, the two surveys are closely coordinated to provide 
an optimized sampling over wavelength. 

In Gruppioni et al. (2010), we started to determine the evolution 
with redshift of the galaxy and AGN LF in the far-IR domain by 
exploiting the PEP data obtained in GOODS-N by the PEP Science 
Demonstration Programme (SDP). Here, we extend the analysis to 
the wider and shallower fields - COSMOS and ECDFS - and to 
the deepest field - GOODS-S - observed by PEP, and we also take 
advantage of the HerMES sub-mm data in the same fields to derive 
improved SED classifications and accurate total IR luminosities for 
our sources. This allows us to have statistically significant samples 
of IR galaxies at different redshifts and over a broad range of lumi- 
nosities, to make a detailed study of the LF at several z intervals, all 
the way from z = 0 to z — 4. The measure of the total IR luminosity 
obtained by integrating the SEDs, well constrained over the entire 
mid-IR and far-IR domain (and also in the sub-mm thanks to the 
available SPIRE data), allows us to derive the total IR LF and its 
evolution directly from far-IR data for unbiased samples selected 
at wavelengths close to the peak of dust emission. Moreover, the 
availability of deep multiwavelength catalogues in the PEP fields 
is crucial for analysing the SEDs, obtaining ^-corrections and total 
IR luminosities, and classifying the PEP sources into different IR 
populations, in order to separately study their LFs and evolutionary 
behaviour. This is the first study ever based on such a statistically 
wide and deep far-IR sample, to be able to provide LFs for different 
IR populations of galaxies and AGN. Here, the evolution of the 
far-IR and total IR LFs (and luminosity density, hereafter pi R ) is 
derived up to unprecedented high redshifts (~4) both globally (e.g. 
for all the populations together) and separately, for each SED class. 

Despite the abundance of information available in the literature 
about the stellar mass function (MF; Fontana et al. 2004; Pozzetti 
et al. 2010; Ilbert et al. 2010; Dominguez-Sanchez et al. 201 1), very 
little is known about the corresponding total IR LFs and SF rate 
(SFR) densities at different masses (an attempt based on Spitzer 
data was made by Perez-Gonzalez et al. 2005). From stellar MF 
studies, one finds a clear increase with z of the relative fraction 
of massive [log(M/MQ >11] star-forming objects, starting to con- 
tribute significantly to the massive-end of the MFs at z > 1 (Fontana 
et al. 2004; Ilbert et al. 2010). Their evolution and contribution to 
the total SFR history are however still uncertain, since only a few 
studies have tried to reconstruct the evolution of the SF history 
of massive objects from optical/near-IR or mid-IR surveys (Juneau 
et al. 2005; Perez-Gonzalez et al. 2005; Santini et al. 2009; Fontanot 
et al. 2012) but none from far-IR- selected surveys (providing a more 
direct indicator of the galaxy SF activity). In this work, we have de- 
rived the IR LF and density in three different mass ranges [from 
log(M/MQ) = 8.5 to log(M/MQ) = 12], extending previous studies 
(limited to z = 1.8-2 for the most massive galaxies) to z ~ 4. 

Finally, our PEP data sets have allowed us to quantify the relative 
contribution of the two main modes of SF (a relatively steady one 
in disc-like galaxies, defining a tight SFR-stellar mass sequence, 
and a starburst mode in outliers) to the total IR LF and pi R in three 
redshift intervals (0.8 < z < 1.25, 1.25 < z < 1.8 and 1.8 < z < 
2.2) and to test the SED-classes belonging to each mode. 

This paper is structured as follows. The PEP Survey with the 
far-IR and multiwavelength data, together with the SED character- 
ization and redshift distribution of the PEP sources, are described 
in Section 2. The LFs (rest-frame 35-, 60- , 90- qm and total IR) and 


their evolution [derived for different SED-classes, mass and spe- 
cific SFR (sSFR) intervals] are discussed in Section 3. In Section 4, 
we present the number and IR luminosity densities of the different 
galaxy types, while in Section 5, we discuss our results. In Section 
6, we present our conclusions. 

Throughout this paper, we use a Chabrier initial mass func- 
tion (IMF) and we assume a A cold dark matter cosmology with 
H 0 = 71 km s -1 Mpc -1 , Q m = 0.27 and £2 A = 0.73. 

2 THE DATA 

The PEP fields where we computed the LFs are COSMOS, 2 deg 2 
observed down to 3cr depths of ~5 and 10.2 mJy at 100 and 160 qm, 
respectively; ECDFS, ~700 arcmin 2 down to 3 a depths of ~4.5 and 
8.5 mJy at 100 and 160 qm, respectively; GOODS-N, ~300 arcmin 2 
to 3 and 5.7 mJy at 100 and 160 qm, respectively; and GOODS- 
S, ~300 arcmin 2 to 1.2, 1.2 and 2.4 mJy at 70, 100 and 160 qm, 
respectively. Our reference samples are the blind catalogues at 70 
(in GOODS-S only), 100 and 160 qm to the 3<j level, which con- 
tain 373 (all in GOODS-S), 7176 (GOODS-S: 717, GOODS-N: 
291, ECDFS: 813, COSMOS: 5355) and 7376 (GOODS-S: 867, 
GOODS-N: 316, ECDFS: 688, COSMOS: 5105) sources at 70, 100 
and 160 qm, respectively. We refer to Berta et al. (2010, 2011) for 
a detailed description of the data catalogues and source counts. 

2.1 Multiwavelength identification 

The PEP fields benefit from an extensive multi wavelength coverage. 
We have therefore associated our sources to the ancillary catalogues 
by means of a multiband likelihood ratio technique (Sutherland 
& Saunders 1992; Ciliegi et al. 2001), starting from the longest 
available wavelength (160 qm, PACS) and progressively match- 
ing 100 qm (PACS), 70 qm (PACS, GOODS-S only) and 24 qm 
(Spter/Multi-band Imaging Photometer - MIPS). In the GOODS- 
S field, we have associated with our PEP sources the 24- qm cat- 
alogue by Magnelli et al. (2009) that we have matched with the 
optical+near-IR+IRAC MUltiwavelength Southern Infrared Cat- 
alogue (MUSIC) catalogue of Grazian et al. (2006), revised by 
Santini et al. (2009), which includes spectroscopic and photometric 
redshifts. To maximize the fraction of identifications, we limited our 
study to the area covered by the MUSIC catalogue (~196 arcmin 2 ), 
obtaining 233, 468 and 492 sources at 70, 100 and 160 qm, re- 
spectively, with flux density greater than the flux limits reported 
above (all with either spectroscopic or photometric redshifts). In 
the GOODS-N field, as described in Berta et al. (2010, 2011) and 
Gruppioni et al. (2010), a PSF-matched multiwavelength catalogue 1 
was created, including photometry from the far-UV ( GALEX) to the 
mid-IR (, Spitzer ). As in GOODS-S, to maximize the identifications, 
we limited our study in GOODS-N to the area covered by the Ad- 
vanced Camera for Surveys (ACS) on the Hubble Space Teleacope 
(~150 arcmin 2 ), obtaining 176 and 186 sources with flux density 
greater than the flux limit at 100 and 160 qm, respectively (all with 
redshifts). We have matched our sources in the ECDFS with the 
Multi wavelength Survey by Yale-Chile by Cardamone et al. (2010), 
obtaining 687 sources at 100 qm and 625 sources at 160 qm (578 
and 547 with redshifts, ~45 per cent spectroscopic). Finally, in 
COSMOS, we have matched our catalogue with the deep 24- qm 
sample of Le Floc’h et al. (2009) and with the IRAC-based catalogue 
of Ilbert et al. (2010), including optical and near-IR photometry and 

1 Publicly available at http://www.mpe.mpg.de/ir/Research/PEP/public_ 
data_releases.php. 
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photometric redshifts. After the removal of PEP sources within 
flagged areas of the optical and/or IRAC COSMOS catalogues, we 
ended up with two catalogues consisting of 41 10 and 4118 sources, 
with flux densities >5.0 and >10.2 mJy at 100 and 160 jam, re- 
spectively (3817 and 3849 with either spectroscopic or photometric 
redshifts). Throughout this paper and specifically for the SED fits 
described in Section 2.2, we adopt these spectroscopic or rest-frame 
UV to near-IR photometric redshifts for the various fields. 

The HerMES extragalactic survey (Oliver et al. 2012) performed 
coordinated observations with SPIRE at 250, 350 and 500 qm in 
the same fields covered by PEP. In particular, in HerMES a prior 
source extraction was performed using the method presented in 
Roseboom et al. (2010), based on MIPS-24 jam positions. The 24- 
|am sources used as priors for SPIRE source extraction are the same 
as those associated with our PEP sources through the likelihood 
ratio technique. We have therefore associated the HerMES sources 
with the PEP sources by means of the 24- jam sources matched 
to both samples. For most of our PEP sources (~87 per cent on 
average, though the SPIRE identification rates vary from ~84 per 
cent in GOODS-S to 96 per cent in ECDFS), we found a >3cr 
SPIRE counterpart in the HerMES catalogues. 

2.2 Galaxy classification 

We made use of all the available multiwavelength data to derive the 
SEDs of our PEP sources, which we interpreted and classified by 
performing a x 2 fit (using the le phare code; 2 Arnouts et al. 2002 
and Ilbert et al. 2006) with the semi-empirical template library of 
Polletta et al. (2007), representative of different classes of IR galax- 
ies and AGN. To this library we added some templates modified in 
their far-IR part to better reproduce the observed Herschel data (see 
Gruppioni et al. 2010), and three starburst templates from Rieke 
et al. (2009). If required to improve the fit, different extinction val- 
ues [E(B — V) from 0.0 to 0.5] have been applied to the templates, 
by letting the code free to choose the most suitable extinction curve. 
The considered set of templates included SEDs of elliptical galax- 
ies of different ages, lenticular, spirals (from Sa to Sdm), starburst 
galaxies (SB), type 1 QSOs, type 2 QSOs, Seyferts, LINERs and 
composite Ultra Luminous Infrared Galaxies (ULIRGs, containing 
both starburst and obscured AGN component), in the wavelength 
range 0.1-1000 qm. The latter templates are empirical ones created 
to reproduce the SEDs of the heavily obscured AGN. Two of these 
SEDs [the broad absorption-line QSO Markarian 231 (Berta 2005) 
and the Seyfert 2 galaxy IRAS 19254—7245 South (Berta et al. 
2003)] are similar in shape, containing a powerful starburst com- 
ponent, mainly responsible for their far-IR emission, and an AGN 
component that contributes to - and dominates - the mid-IR (Farrah 
et al. 2003), reproducing the SEDs of ‘obscured’ AGN regardless 
of their optical spectra (i.e. broad or narrow lines in the optical; 
Gruppioni et al. 2008). Hereafter, we will refer to this class of 
templates and to the sources reproduced by them as to type 2 AGN 
(AGN2). Three other empirical templates, reproducing the observed 
SEDs of nearby ULIRGs containing an obscured AGN (i.e. IRAS 
20551-4250; IRAS 22491-1808; NGC 6240) have been associ- 
ated with the Seyfert 1.8/2, LINER ones, since they all contain 
an AGN, but this AGN does not dominate the observed energetic 
output at any wavelength (from UV to far-IR/sub-mm), showing 
up just in the range where the host galaxy SED has a minimum 
(i.e. the mid-IR). The AGN in these objects is either obscured or 
of low luminosity. We refer to this class as to star-forming galaxies 

2 Available at http://www.cfht.hawaii.edu/~arnouts/LEPHARE/lephare.html. 


containing an AGN (SF-AGN), since their IR luminosity is largely 
dominated by SF. 

In our analysis, we make the basic assumption that the SED 
shapes seen at low redshifts are also able to represent the higher 
redshift objects. In any case, to further increase the range of SEDs 
in the fit, we have applied additional extinction with different ex- 
tinction curves to our templates. All SED fits adopt spectroscopic 
or photometric redshifts described in Section 2.1. 

The template library used to fit our data contains a finite number 
of SEDs (38), representative of given classes of local infrared ob- 
jects, which do not vary with continuity from one class to another 
(there are large gaps in the parameter space). Therefore, the quality 
of the fit depends not only on the photometric errors, but also on 
the template SED uncertainties. For this reason, in our fitting proce- 
dure, in addition to the photometric errors on data, we need to take 
into account also the uncertainties due to the template SEDs dis- 
cretization and additional extinction. To do this, we have proceeded 
as described in detail by Gruppioni et al. (2008) and summarized as 
follows. First, we have run le phare on our PEP SEDs considering 
the nominal errors from catalogues, computing the distributions of 
the (Subject Ne m pi ate) band/ (^ object )band values in each of the consid- 

ered photometric band (where Subject and object are the flux density 
and the relative error of the source, and template the flux density 
of the template in the considered band), iteratively increasing the 
photometric errors until we have obtained a Gaussian distribution 
with a ~ 1. This corresponds to reduced yf distributions peaked 
around 1 (as expected in the case of good fit). With the new photo- 
metric uncertainties (on average, significantly increased mainly in 
the optical/near-IR and SPIRE bands), we have run le phare on our 
sources for the second time, obtaining what we have taken as the 
final SED-fitting results. 

The majority of our PEP sources are reproduced by templates 
of normal spiral galaxies (spiral), SB galaxies (starburst) and 
Seyfert2/1 ,8/LINERS/ULIRGs-f AGN (SF-AGN), although differ- 
ent classes prevail at different redshifts and luminosities. The spiral 
SEDs show no clear signs of enhanced SF or nuclear activity (see 
Fig. 1), the far-IR bump being characterized by relatively cold dust 
(Tdust ~ 20 K). On the other hand, SB templates are characterized 
by warmer (r dust ~ 40-45 K), more pronounced far-IR bumps and 
significant UV extinction, indicative of intense SF activity. Tem- 
plates of star-forming galaxies containing either a low-luminosity 
or obscured AGN (SF-AGN) are characterized by a ‘flattening’ in 
the 3-10 qm spectrum (suggesting detection of an AGN in the 
wavelength range where the host galaxy SED has a minimum) and 
a far-IR bump dominated by SF, which is intermediate (in terms of 
both energy and r dust ) between spirals and SBs. Although they can 
be considered as star-forming galaxies at the wavelengths relevant 
to this work, we prefer to refer to them as SF-AGN throughout the 
paper, to keep in mind that they probably contain an AGN, whose 
presence, though not dominant in the far-IR, might be very impor- 
tant for analysis in other bands (e.g., in the X-rays or the mid-IR). 

We note that the well-studied high-z (—2.3) SED of the strongly 
lensed sub-mm galaxy SMM J2135— 0102, known as ‘the Cosmic 
Eyelash’ (Ivison et al. 2010; Swinbank et al. 2010), best fitted with 
our procedure by an extinct [E{B — V) ~ 0.2] IRAS 22491 — 1808 
template (though rather poorly in the near-IR), does not represent 
the bulk of our population at high-z (>1.5), whose SEDs are indeed 
well reproduced by our library of templates. 

In fact, the considered template set provides very good fits to 
the SEDs of our PEP sources. In Fig. 1, we show the rest-frame 
SEDs (black dots) of the PEP sources belonging to the different 
‘broad’ SED classes (spiral, starburst, SF-AGN, AGN1 and AGN2), 
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Figure 1 . Observed rest-frame SEDs of the PEP sources (black dots) divided by population (as shown in the plot) and normalized to the K s band. The more 
representative templates for each SED-class have been overplotted in different colours. 


compared to the template SEDs of those classes normalized to the 
^ s -band flux density. In Table 1, we report the fraction of sources 
belonging to each SED class: we find that in all the fields ~41(38) 
per cent of the 100(1 60)- jam sources are reproduced by a spiral 
template SED, 7(7) per cent with a starburst template SED, 45(48) 
per cent with an SF-AGN template SED, 2(3) per cent with an AGN2 
SED and 5(4) per cent with an AGN 1 SED. We note that the fraction 
of SF-AGN derived in this work is in agreement with results from 
mid-IR spectroscopy (with Spitzer - InfraRed Spectrograph) of local 
star-forming galaxies from the Spitzer Infrared Nearby Galaxies 
Survey (SINGS) sample by Smith et al. (2007), who found that 
~50 per cent of local galaxies (though of lower luminosities than 
ours) do harbour low-luminosity AGN (of LINER or Seyfert types). 
Recently, Sajina et al. (2012) found an even higher fraction (~70 per 
cent) of objects hosting an AGN in the mid-IR (, Spitzer 24- jam) 
selected samples (~23 per cent AGN-dominated and ~47 per cent 
showing both AGN and starburst activity). However, since the far-IR 
SED of the SF-AGN is dominated by SF and at these wavelengths 


resembles either starburst or spiral galaxy templates, we have also 
divided the SF-AGN class into SF-AGN(SB) and SF-AGN(Spiral) 
sub-classes, based of their far-IR/near-IR colours (e.g. Sm/Si^) 
and SED resemblance (apart from the rest-frame mid-IR flattening, 
which is detected in all of the SF-AGN SEDs). Specifically, galaxies 
best- fitted by the Seyfert 2/Seyfert 1.8 templates (either the original 
ones from Polletta et al. 2007 or those modified by Gruppioni et al. 
2010) have been classified as SF-AGN(Spiral), while galaxies best- 
fitted by the NGC 6240, IRAS 205511-4250 or IRAS 22491-1808 
templates have been classified as SF-AGN(SB). The number of 
sources belonging to the former and the latter sub-classes are also 
reported in Table 1 as additional information. 

2.3 Redshift distribution 

A large number of spectroscopic redshifts have been measured in 
the GOODS, ECDFS and COSMOS regions. In the GOODS-S and 
ECDFS area, a collection of more than 5000 spectroscopic redshifts 
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Table 1 . SED classification. 


Field 

Spiral 

Starburst 

SF-AGN 

AGN2 

AGN1 

SF-AGN(SB) 

SF-AGN(Spiral) 

GOODS-S 

N 

N 

N 

N 

N 

N 

N 

70 pirn 

53 

22 

142 

5 

12 

26 

116 

100 pm 

(23 per cent) 
117 

(9 per cent) 
60 

(61 per cent) 
250 

(2 per cent) 
10 

(5 per cent) 
31 

54 

96 

160 pm 

(25 per cent) 
123 

(13 per cent) 
55 

(53 per cent) 
277 

(2 per cent) 
11 

(7 per cent) 
26 

73 

204 


(25 per cent) 

(11 per cent) 

(56 per cent) 

(2 per cent) 

(6 per cent) 



GOODS-N 

N 

N 

N 

N 

N 

N 

N 

100 pm 

68 

20 

78 

7 

3 

21 

57 

160 pm 

(39 per cent) 
67 

(11 per cent) 
21 

(44 per cent) 
85 

(4 per cent) 
10 

(2 per cent) 
3 

21 

64 


(36 per cent) 

(11 per cent) 

(46 per cent) 

(5 per cent) 

(2 per cent) 



ECDFS 

N 

N 

N 

N 

N 

N 

N 

100 pm 

253 

49 

245 

8 

23 

83 

162 

160 pm 

(44 per cent) 
233 

(9 per cent) 
49 

(42 per cent) 
231 

(1 per cent) 
12 

(4 per cent) 
22 

99 

132 


(43 per cent) 

(9 per cent) 

(42 per cent) 

(2 per cent) 

(4 per cent) 



COSMOS 

N 

N 

N 

N 

N 

N 

N 

100 pm 

1637 

232 

1689 

76 

183 

580 

1109 

160 pm 

(43 per cent) 
1483 

(6 per cent) 
243 

(44 per cent) 
1847 

(2 per cent) 
103 

(5 per cent) 
173 

111 

1070 


(39 per cent) 

(6 per cent) 

(48 per cent) 

(3 per cent) 

(4 per cent) 



TOTAL 

N 

N 

N 

N 

N 

N 

N 

100 pm 

2075 

361 

2262 

101 

240 

738 

1424 

160 pm 

(41 per cent) 
1906 

(7 per cent) 
368 

(45 per cent) 
2440 

(2 per cent) 
136 

(5 per cent) 
224 

970 

1470 


(38 per cent) 

(7 per cent) 

(48 per cent) 

(3 per cent) 

(4 per cent) 




are available (Cristiani et al. 2000; Croom, Warren & Glazebrook 
2001; Bunker et al. 2003; Dickinson et al. 2004; Stanway et al. 
2004; Strolger et al. 2004; Szokoly et al. 2004; van der Wei et al. 
2004; Doherty et al. 2005; Le Fevre et al. 2005; Mignoli et al. 
2005; Vanzella et al. 2008; Popesso et al. 2009; Santini et al. 2009; 
Balestra et al. 2010; Cooper et al. 2012). In the GOODS-N area, 
more than 2000 spectroscopic redshifts come from various obser- 
vations (Cohen et al. 2000; Cowie et al. 2004; Wirth et al. 2004; 
Barger, Cowie & Wang 2008). Finally, in COSMOS we could use 
a collection of ~3000 spectroscopic redshifts from either the pub- 
lic zCOSMOS bright data base or the non-public zCOSMOS deep 
data base (Lilly et al. 2007, 2009). For the PEP sources without 
spectroscopic redshift available, we have adopted the photometric 
redshifts derived from multi- wavelength (UV to near-IR) photom- 
etry by different authors in the different fields, as mentioned in 
Section 2.1. In the GOODS-S field, the MUSIC photometric red- 
shift catalogue (Grazian et al. 2006; Santini et al. 2009) provided 
photo-z s for most of our PEP sources without spectroscopic data, 
while in the GOODS-N field, photo-zs were obtained by Berta et al. 
(2010) for almost all the PEP sources within the ACS area. The 
Cardamone et al. (2010) and Ilbert et al. (2009) catalogues provided 
photometric redshifts for a large fraction of the PEP sources in the 
ECDFS and COSMOS areas, respectively. When considering both 
the spectroscopic and photometric redshifts, in our PEP fields the 
redshift incompleteness is very low. In particular, in the GOODS-S 
field we have either a spec-z or a photo-z for ~100 per cent of the 
PEP sample within the MUSIC areas (~80 per cent spectroscopic, 
though most of them lie at z < 2.5; see Berta et al. 2011). In the 


GOODS-N field, we have a redshift completeness of ~100 per cent 
of sources (70 per cent spectroscopic) within the ACS area. In the 
ECDFS and COSMOS fields, we have a redshift completeness of 
88 and 93 per cent, respectively (45 and 40 per cent spectroscopic). 

The uncertainty in the photometric redshifts has been evaluated 
by means of a comparison with the available spec-z s by the dif- 
ferent authors providing photo-z catalogues in the PEP fields. In 
particular, Berta et al. (2011) have compared the photometric and 
the available spectroscopic redshifts in GOODS-S, GOODS-N and 
COSMOS, finding a fraction of outliers, defined as objects having 
Az/(1 + Zspec) > 0.2, of ~2 per cent for sources with a PACS 
detection. Most of these outliers are sources with few photomet- 
ric points available, or SEDs not well reproduced by the available 
templates. The median absolute deviation of the Az/(1 + z spe c) 
distribution in the three fields analysed by Berta et al. (2011) is 
0.04 for the whole catalogue, and 0.038 for PACS-detected objects. 
In GOODS-S, Grazian et al. (2006) found excellent agreement be- 
tween photometric and spectroscopic redshifts over the fully ac- 
cessible redshift range 0 < z < 6 (<r[Az/(l + z)] — 0.045), with 
a very limited number of catastrophic errors. In COSMOS, Ilbert 
et al. (2010) estimated the photometric redshift uncertainties of 
their 3.6- jam catalogue matched with the COSMOS photo-z mul- 
tiwavelength catalogue of Ilbert et al. (2009), finding <r[Az/(l + 
z)] = 0.008 (and < 1 per cent of catastrophic failures) at i^ B < 22.5, 
<r[Az/(l + z)] = 0.011 at 22.5 < /+ B < 24 and <r[Az/(l + z)] = 
0.053 at 24 < i^ B < 25. In the ECDFS, by comparing non-X-ray 
sources with high-quality spectroscopic redshifts, Cardamone et al. 
2010 found cr[Az/(l + z)] = 0.008 to z ~ 1.2 = 0.027 at 1.2 < 
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z <3.7 and = 0.016 at z > 3.7. Note that we have checked all the 
z > 2.5 photometric redshifts through le phare, assigning the le 
PHARE-derived value in case of significant disagreement with that 
from the catalogue (though most resulted in very good agreement). 
The fractions of spectroscopic redshifts in the 2.5 < z < 3.0 interval 
amount to just ~6 and ~25 per cent in COSMOS and GOODSS, 
respectively, dropping to ~4 and ~6 per cent at 3.0 < z < 4.2. We 
note that from our comparison between photo- and spec-z (when 
available) we find a general good agreement in all fields. 

Photometric redshift errors may, in principle, affect the shape of 
the LF at the bright end: by scattering objects to higher redshifts, 
they make the steep fall-off at high luminosities appear shallower 
(e.g. Drory et al. 2003). To study the impact of redshift uncertain- 
ties on the inferred infrared LF, we have performed Monte Carlo 
simulations, as discussed in detail in Section 3.3. It is indeed very 
difficult to estimate the effect of catastrophic failures at z > 2.5, 
where mainly photometric redshifts are available and no reliable 
spectroscopic data can be used to validate them. Moreover, for the 
limited high-z samples with spectroscopic information, different 
results are found in the different fields: i.e. in GOODS-S, Berta 
et al. (2011) found that up to 25 per cent of PACS detected sources 
above z ~ 2 might have a higher than real photometric redshift, 
while in COSMOS the catastrophic failures (20 per cent) found by 
Ilbert et al. (2009) for MIPS-selected sources at 2 < z < 3 were 
mostly for photo-z smaller than spectroscopic ones. In addition to 
that, sometimes high-z spectroscopic redshifts can be even more 
uncertain than photometric ones and great care must be taken when 
selecting spec-zs for comparison (i.e. we need to choose those with 
high-quality flags). For all these reasons, we limited our analysis 
of photo-z uncertainties to the Monte Carlo simulations described 
in Section 3.3, without trying to derive uncertainties also due to 
catastrophic failures. In Section 3.2 we also note that the different 
(far-IR) photo-z approach of Lapi et al. (2011) produces consistent 
LF results in the common part of the parameter space. 

The median redshift of the 70- pm sample in GOODS-S is 
£med(70) = 0.67 (the mean is (zb o = 0.86), while those of the 
100- and 160- pm samples are different in each of the fields, given 
the different flux density depths reached by PEP in each area. In 
Table 2, we report the median and the mean redshifts found for 
the different fields (and for the combined sample) at the different 
selection wavelengths. As expected, GOODS-S reaches the highest 
redshifts, while the surveys in the COSMOS and ECDFS fields are 
shallower and sample lower redshifts. On average, the 160- pm se- 
lection favours higher redshifts than the 100- pm one (see also Berta 
et al. 2011). 

In Fig. 2, we show the redshift distribution of the PEP sources 
selected at 100 and 160 pm in the four different fields. The black 
solid histogram is the total redshift distribution in the field (one 
for each row of the plot), while the filled histograms in different 
colours represent the redshift distributions of the different popula- 
tions (green, spiral; cyan, starburst; red, SF-AGN; magenta, AGN2; 


Table 2. Average redshifts in the PEP Survey fields. 


Field 

70 pm 

100 pm 

160 pm 


Zmed 

(z) 

Zmed 

(z) 

Zmed 

(z) 

GOODS-S 

0.67 

0.86 

0.85 

1.07 

0.98 

1.16 

GOODS-N 



0.73 

0.84 

0.84 

0.94 

ECDFS 



0.66 

0.76 

0.69 

0.85 

COSMOS 



0.59 

0.74 

0.70 

0.88 

ALL FIELDS 

0.67 

0.86 

0.64 

0.78 

0.73 

0.91 


blue, AGN1). The line-filled dashed histograms shown in the spi- 
ral and starburst panels represent the redshift distributions of the 
SF-AGN(Spiral) and SF-AGN(SB) sub-classes, respectively. In ad- 
dition to the principal redshift peak, in GOODS-S a secondary peak 
centred at z ~ 2 is clearly visible at both 100 and 160 pm. A similar 
result has been shown and discussed also by Berta et al. (2011), 
while an extensive analysis of PACS GOODS-S large-scale struc- 
ture at z = 2-3 and of a z ~ 2.2 filamentary overdensity have been 
presented by Magliocchetti et al. (2011). 

3 THE LUMINOSITY FUNCTION 

The sizes and depths of the PEP samples are such as to allow a di- 
rect and accurate determination of the far-IR LF in several redshift 
bins, from z — 0 up to z ~ 4. PEP+HerMES is the unique Herschel 
survey to allow such analysis over such a wide redshift and lumi- 
nosity range, sampling both the faint and bright ends of the far-IR 
and total IR LFs with sufficient statistics. Because of the redshift 
range covered by PEP, we would need to make significant extrap- 
olations in wavelength when computing the rest-frame LFs at any 
chosen wavelengths. In order to apply the smallest extrapolations 
for the majority of our sources, we choose to derive the far-IR LFs 
at the rest-frame wavelengths corresponding to the median redshift 
of each sample. Given the median redshift of the 70- pm sample 
in GOODS-S (~0.67, see Table 2), we use that sample to derive 
the rest-frame LF at 35 pm. With the 100- and 160- pm PEP sam- 
ples (whose median redshifts are ~0.64 and 0.73, respectively), we 
derived the rest-frame LFs at 60 and 90 pm. Note that, given the ex- 
cellent multiwavelength coverage available for most of our sources 
(thanks also to the HerMES data available in all the PEP fields and 
providing reliable counterparts for most of our PEP sources), their 
SEDs are very well determined from the UV to the sub-mm. The 
extrapolations are therefore well constrained by accurately defined 
SEDs, even at high redshifts (i.e. at z ~ 3.5 the rest-frame 90-pm 
luminosity corresponds to A observed ~ 400 pm, which is still in the 
range covered by HerMES). 

3.1 Method 

The LFs are derived using the 1 / V max method (Schmidt 1968). This 
method is non-parametric and does not require any assumptions on 
the LF shape, but derives the LF directly from the data. We have 
first derived the LFs in each field separately, in order to check for 
consistency and to test the role of cosmic variance. Successively, we 
have made use of the whole data sets to derive the monochromatic 
and total IR LFs, by means of the Avni & Bahcall (1980) method 
for coherent analysis of independent data sets. We have divided 
the samples into different redshift bins, over the range 0 < z < 4, 
selected to be almost equally populated, at least up to z ~ 2.5. In 
each redshift bin, we have computed the comoving volume available 
to each source belonging to that bin, defined as V max = V Zmax — V Zmin , 
where z ma x is the minimum between the maximum redshift at which 
a source would still be included in the sample given the limiting 
flux of the survey (different for each field) and the upper boundary 
of the considered redshift bin, while z m in is just the lower boundary 
of the considered redshift bin. 

When combining the four samples, we have constructed 
a complete sample over the whole GOODS-S+GOODS-N+ 
ECDFS(— GOODS-S)+COSMOS region, including all the ob- 
served objects (see details in Section 3.2). The depth of the sample 
is not constant throughout the region, but an object with a given flux 
density (included in the list of observed objects irrespective of the 
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Figure 2. Redshift distributions of 100- pm (top) and 160- pm (bottom) sources in the four PEP fields (first row from top, GOODS-S; second row, GOODS-N; 
third row, ECDFS; bottom row, COSMOS) to different limiting fluxes (as shown in the plot). The redshift distributions of the five different populations have 
been plotted in different colours (green, spiral; cyan, starburst; red, SF-AGN; magenta, AGN2; blue, AGN1) and compared to the total distribution (black 
solid histogram). The line-filled dashed histograms shown in the spiral and starburst panels represent the redshift distributions of the SF-AGN(Spiral) and 
SF-AGN(SB) sub-classes, respectively. 
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identity of its parent sample) can a priori be found in one (or more) 
region if its redshift is <zJ^(Siimit) of that region (e.g. sources de- 
tected in the COSMOS area are detectable over the whole joint area, 
while the fainter sources detected in GOODS -S are detectable in 
GOODS-S only). The maximum volume of space which is available 
to such an object to be included in the joint sample is then defined 
by 
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where V^ ax (with fid = GS, GN, E, C corresponding to GOODS-S, 
GOODS -N, ECDFS and COSMOS, respectively) is the comoving 
volume available to each source in that field, in a given redshift bin, 
while is the solid angle subtended by that field sample on the 
sky. 

For each luminosity and redshift bin, the LF is given by 




1 


w i x V m 


( 2 ) 


where V max , i is the comoving volume over which the ith galaxy 
could be observed, A L is the size of the luminosity bin and wi is 
the completeness correction factor of the ith galaxy. These com- 
pleteness correction factors are a combination of the completeness 
corrections given by Berta et al. (2010, 2011), derived as described 
in Lutz et al. (2011), to be applied to each source as a function of its 
flux density, together with a correction for redshift incompleteness 
(for the ECDFS and COSMOS only, see Section 2.3). Since, as 
mentioned in Section 2.3, the redshift incompleteness in the COS- 
MOS and ECDFS areas is independent on PACS flux density, in 
these fields we have applied the corrections regardless of the source 
luminosity and redshift (i.e. by multiplying <3>(L, z) by 1.07 and 
1.14 in COSMOS and ECDFS, respectively). However, the redshift 
incompleteness does not affect our conclusions, since >95 per cent 
of all our sources do have a redshift. 

Uncertainties in the infrared LF values depend on photometric 
redshift uncertainties. To quantify the effects of the uncertainties in 
photometric redshifts on our luminosity functions, we performed a 
set of Monte Carlo simulations, as described in Section 3.3. 


3.2 The rest-frame 35-, 60- and 90- pm LF 

By following the method described above, we have derived the 
35-, 60- and 90- pm rest-frame LFs from the 70- (in GOODS-S 
only), 100- and 160- pm samples, respectively. In order to check 
the consistency between the catalogues and the effects of cosmic 
variance, we have first derived the monochromatic LFs in each field 
separately. Note that, since the 70- pan data are available in the 
GOODS-S field only, to have a better sampling of the LF especially 


at the bright-end, we have also computed the rest-frame 35- pm LF 
from the 100- pm samples in the four fields and compared them 
with that obtained from the 70-pm sample (see Fig 3; Table 3). 
The agreement between the two derivations is very good, implying 
correct extrapolations in wavelength due to the good and complete 
SED coverage. 

We have divided the samples into seven redshift bins: 0.0 < z < 
0.4; 0.4 < z < 0.8; 0.8 < z < L2; 1.2 < z < 1.8; 1.8 < z < 2.5; 2.5 < 
z < 3.5 and 3.5 < z <4.5. The results of the computation of our 35- 
(reported in Table 3), 60- and 90- pm LFs are shown in Figs 3, 4 
and 5, respectively. The LFs in the four different fields appear to be 
consistent with each other within the error bars (± lcr ) in most of the 
common luminosity bins. The COSMOS and GOODS-S Surveys 
are complementary, with the faint end of the LFs being mostly 
determined by data in GOODS-S, and the bright end by COSMOS 
data. After having checked the field-to-field consistency, we have 
combined the 100- and 160- pm samples in all fields, obtaining the 
global rest-frame 60- and 90- pm LFs (reported in Tables 4 and 5, 
respectively). 

The data from each field in each z-bin have been plotted (and 
considered in the combination) only in the luminosity bins where we 
expect our sample to be complete, given that at fainter luminosities 
not all galaxy types are observable (depending on their SEDs; Ilbert 
et al. 2004). For comparison, we overplot the LFs at 35 pm from 
Magnelli et al. (2009, 201 1), the local LFs at 60 pm from Saunders 
et al. (1990) and those at 90 pm from Serjeant et al. (2004) and 
Sedgwick et al. (2011) and at 100 pm from Lapi et al. (2011), 
respectively. The comparison between the 35- pm LF, derived from 
the 70- pm PEP sample in GOODS-S, and the results of Magnelli 
et al. (2009, 201 1), based on a 24- pm prior extraction and stacking 
analysis on Spitzer maps, shows very good agreement, both with 
the data and with the double power-law fit. The 1.8 < z < 2.5 PEP 
LF is consistent within ±1 o with the Magnelli et al. (2011) data 
points, though the power-law fit at bright L 35 (> 10 12 Lq) seems to 
be slightly lower than our data. At z > 2.5 no comparison data from 
Spitzer are available, while our LF derivation can provide hints 
of evolution at the bright end of the LF. In the common redshift 
intervals (between z ~ 1 and 3.5), our 90- pm LF is in very good 
agreement with the 1 00- pm Lapi et al. (20 11) derivation from the H- 
ATLAS survey (although their redshift bins are somewhat different 
than ours: 1.2 < z < 1.6; 1.6 < z < 2.0; 2.0 < z < 2.4 and 2.4 < 
z < 4.0) and with the previous PEP-SDP derivation (Gruppioni 
et al. 2010). The consistency between our 90- pm LFs and the Lapi 
et al. (201 1) ones (derived from a different sample, using a different 
template SED to fit the data and a rest-frame far-IR based method 
to derive photometric redshifts) gives us confidence that, at least 
up to z ~ 3.5, our computation is not significantly affected by 
incompleteness or by photo-z uncertainties. 


3.3 The total infrared LF 

We integrate the best-fitting SED of each source over 8 < A. rest < 
1000 pm to derive the total IR luminosities (Li R = L[8-1000 pm]) 
in 11 redshift bins (0.0-0.3; 0.3-0.45; 0.45-0.6; 0.6-0.8; 0.8-1.0; 
1.0-1. 2; 1.2-1. 7; 1.7-2.0; 2.0-2.5; 2.5-3.0 and 3.0-4.2), selected 
to be almost equally populated, at least up to z ~ 2.5. Our approach 
is similar to that of other studies based on mid-IR- selected galaxy 
samples (e.g. Le Floc’h et al. 2005; Rodighiero et al. 2010a; Mag- 
nelli et al. 2011), but this is the first time that the SEDs have been 
accurately constrained by sufficiently deep data in the far-IR/sub- 
mm domain and not simply extrapolated from the mid-IR to the 
longer wavelengths or from average flux density ratios. 
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Figure 3. Rest- frame 35- jam LF estimated with the 1 /V max method from the PEP 70- jam sample in the GOODS-S field (purple open circles) and independently 
in the four PEP fields from the 100- pm-selected samples (red filled circles, GOODS-S; green filled triangles, GOODS-N; orange filled stars, ECDFS; blue 
filled squares, COSMOS). The error bars in the data points represent the Poissonian uncertainties. For comparison, we also plot the determination of Magnelli 
et al. (2011) at 1.2 < z < 2.5, shown as green open squares, and the double power-law fit of Magnelli et al. (2009, 2011), shown as a green dashed line. The 
black dashed line represents the z = 0 determination of Magnelli et al. (2009). 


As mentioned in Section 2.3, we have studied the impact of red- 
shift uncertainties on the total IR LF by performing Monte Carlo 
simulations. As test cases, we used the COSMOS and GOODS-S 
samples (which are basically defining the bright and faint ends of 
the LF in an almost complementary way), and we checked the ef- 
fect on the total IR LF. We iterated the computation of the total IR 
LF by each time varying the photometric redshift of each source 
(assigning a randomly selected value, according to a Gaussian dis- 
tribution, within the 68 per cent confidence interval). Each time, 
we then recomputed the monochromatic and total IR luminosities, 
as well as the V max value, but we did not perform the SED fitting 
again, keeping the previously found best-fitting template for each 
object (the effect on the ^-correction is not relevant in the far-IR 
wavelength range). The results of this Monte Carlo simulation are 
reported in Fig. 6, where we show the total IR LFs derived in each 
of the PEP fields independently: the red and blue filled circles rep- 
resent the estimates of the GOODS-S and COSMOS LFs, with their 
range of values derived with 20 iterations by allowing a change in 
photo-z represented by the pink and sky-blue shaded areas, respec- 
tively. The comparison shows that the effect of the uncertainty of 
the photometric redshifts on the error bars is slightly larger than 
the simple Poissonian value (1/VA), and affects mainly the lower 
and the higher redshift bins (especially at low and high luminosi- 
ties). Using these Monte Carlo simulations, we find no evident 
systematic offsets caused by the photometric redshift uncertainties 
(see Fig. 6). This is due to the very accurate photometric redshifts 
available in these fields. For the total uncertainty in each luminos- 
ity bin in GOODS-S and COSMOS, we have therefore assumed 
the dispersion given by the Monte Carlo simulations (as shown in 


Fig. 7). We note that at the higher z (>2.5), where we must rely on 
a majority of photometric redshifts, the true uncertainties (taking 
into account also catastrophic failures or incompleteness effects) 
might be larger than derived with simulations. The unavailability of 
a ‘true’ comparison sample (i.e. a large comparison sample with ac- 
curate spectroscopic redshifts and fully representative for the PACS 
flux selection) at high z does not allow us to properly quantify this 
statement. 

In Fig. 7, the total IR LFs obtained by combining the 160-qm- 
selected samples with the Avni & Bahcall (1980) technique is plot- 
ted and compared with other derivations available in the literature. 
The total IR LF of Sanders et al. (2003) is plotted as a local refer- 
ence, in addition to the LFs of Le Floc’h et al. (2005), Rodighiero 
et al. (2010a), Caputi et al. (2007) and Magnelli et al. (2009, 2011) 
in various redshift intervals. Globally, data from surveys at differ- 
ent wavelengths agree relatively well over the common z-range. 
No data for comparison are available at z > 2.5, apart from the 
IR LF of sub-mm galaxies from Chapman et al. (2005) and Wall, 
Pope & Scott (2008) at z ~ 2.5, which represent reasonably well 
just the very bright end of the total IR LF. Our derivation is the 
first at such high redshifts, especially in the 3 < z < 4.2 range. 
Note the good agreement between our PEP-based total IR LF and 
the HerMES-based one derived by Vaccari et al. (in preparation), 
shown by the red open squares in Fig. 7. Though consistent within 
the error bars, in the highest redshift bin the HerMES LF is slightly 
higher than ours. Since the 250- qm HerMES selection favours the 
detection of higher redshift sources than the PEP one, it is more 
likely that the PEP LF in the higher-z bin is affected by some 
flux/redshift incompleteness rather than by the presence of low-z 
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Table 3. Rest-frame 35- jam LF. 


k>g(T35/LQ) 

0.0 < z < 0.4 

0.4 < z < 0.8 

0.8 <z < 1.2 

log(0/Mpc 3 dex ! ) 
1.2 <z <1.8 

1.8 <z <2.5 

2.5 <z <3.5 

3.5 <z <4.5 




From the GOODS-S 70- pm sample 




7.9-8. 3 

-1.34 ±0.29 







8. 3-8.7 

-1.62 ±0.04 







8.7-9. 1 

-2.01 ± 0.20 







9. 1-9.5 

-2.38 ±0.14 







9.5-9.9 

-2.63 ±0.12 

-2.26 ± 0.29 






9.9-10.3 

-2.96 ±0.15 

-2.81 ±0.14 






10.3-10.7 

-3.11 ±0.18 

-2.70 ± 0.06 






10.7-11.1 

-3.50 ± 0.43 

-3.23 ± 0.09 

-3.03 ± 0.09 





11.1-11.5 

-3.90 ± 0.43 

-3.98 ± 0.22 

-3.50 ±0.10 

-3.19 ±0.15 




11.5-11.9 


-4.59 ± 0.43 

-4.35 ± 0.25 

-4.14 ±0.17 

-3.82 ± 0.20 



11.9-12.3 




-4.41 ±0.19 

-4.14 ±0.19 



12.3-12.7 





-4.68 ± 0.32 

-4.27 ± 0.30 


12.7-13.1 






-5.27 ± 0.43 

-5.13 ±0.43 

13.1-13.5 







-5.34 ± 0.43 




From the combined 100- pm sample 




1 . 5 - 1. 9 

-1.74 ±0.43 







7.9-8. 3 

-1.64 ±0.23 







8. 3-8.7 

-2.00 ±0.19 







8.7-9. 1 

-2.01 ± 0.09 







9. 1-9.5 

-2.37 ± 0.03 

-2.43 ± 0.37 






9.5-9.9 

-2.53 ± 0.02 

-2.48 ±0.19 






9.9-10.3 

-2.73 ± 0.02 

-2.66 ± 0.08 






10.3-10.7 

-3.14 ±0.03 

-2.56 ± 0.03 

-2.79 ±0.10 





10.7-11.1 

-3.91 ± 0.07 

-3.34 ± 0.02 

-2.91 ± 0.06 





11.1-11.5 

-4.73 ±0.18 

-4.03 ± 0.04 

-3.46 ± 0.03 

-3.13 ±0.07 

-3.24 ± 0.32 



11.5-11.9 

-5.51 ±0.43 

-4.90 ±0.10 

-4.09 ± 0.03 

-3.77 ± 0.04 

-3.38 ±0.10 



11.9-12.3 


-5.60 ± 0.22 

-4.90 ± 0.07 

-4.37 ± 0.03 

-3.98 ± 0.07 

-3.18 ±0.29 


12.3-12.7 



-6.16 ±0.31 

-5.34 ± 0.09 

-4.54 ± 0.04 

-4.30 ±0.12 


12.7-13.1 




-6.46 ±0.31 

-5.48 ± 0.09 

-4.99 ± 0.07 

-4.40 ± 0.39 

13.1-13.5 






-6.03 ±0.15 

-5.91 ± 0.26 

13.5-13.9 






-7.04 ± 0.43 



sources erroneously placed at high-z by incorrect photometric 
redshift assignment (catastrophic failures). The values of our total 
IR LF for each redshift and luminosity bin are reported in Table 6. 


3.4 Evolution 


In order to study the evolution of the total IR LF, we derive a 
parametric estimate of the LF at different redshifts. For the shape 
of the LF, we assume a modified- Schechter function (i.e. Saunders 
et al. 1990), where O(L) is given by 


d>(L)dlogL = O* 



1 

2a 2 


!og?o 



dlogL, 

(3) 


behaving as a power law for L < L* and as a Gaussian in log L for L > 
L*. The adopted LF parametric shape depends on four parameters 
(a, a, L* and <£>*), whose best- fitting values and uncertainties have 
been found using a non-linear least-squares fitting procedure. In 
detail, while in the first z-bin all the parameters have been estimated, 
starting from the second z-bin, the values of a and o have been 
frozen at the values found at lower redshift, leaving only L* and 
free to vary (see Table 7). Note that, in the highest redshift 
bin (3.0 < z < 4.2) we are not able to constrain the LF break, 
while we are up to ~3. Therefore, the results found at this redshift 
are affected by larger uncertainties than the z < 3 ones. However, 
although there is some degeneracy in the values of L* and <£>* at 


3.0 < z < 4.2, the range of allowed value combinations still giving 
a reasonable fit to the three observed data points (constraining the 
bright-end of the LF) is limited and do not significantly affect our 
results. 

In Fig. 8, we plot the total IR LFs at all redshifts with the ±1 o 
Poissonian uncertainty regions (different colours for different z- 
intervals). There is a clear luminosity evolution with redshift, at 
least up to z ~ 3. The apparent ‘fall’ of the z > 3 LF with respect 
to those at lower redshifts, if real, might imply a global negative 
evolution of the IR galaxies and/or AGN starting at z > 3. However, 
as mentioned above, we must point out that in the highest redshift 
bin PACS data could be affected by incompleteness, since with 
increasing redshift (and for intrinsically ‘cold’ sources) the true 
PACS fluxes might fall below the detection limit (i.e. faint sources 
should be missed even if completeness corrections are perfect), 
while, if luminous enough, they can still be detectable by SPIRE. 
This effect is expected to be more relevant in COSMOS, where 
SPIRE data are quite deep relatively to PACS, while it should not 
happen in GOODS-S, where PACS data are very deep compared 
to the SPIRE ones (which are limited by confusion). The high- 
z incompleteness of PACS surveys might also be emphasized by 
the redshift incompleteness of the COSMOS sample, which could 
affect the highest redshift bins more than the lower ones (although 
the redshift incompleteness seems to be independent of redshift; see 
Berta et al. 2011). However, we must point out that a decrease at z > 
2.7-3, similar to that observed in our data, is also observed in the 
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Figure 4. Rest-frame 60- jam LF estimated independently from the 100- jam-selected samples in the four PEP fields (red filled circles, GOODS-S; green filled 
triangles, GOODS-N; orange filled stars, ECDFS; blue filled squares, COSMOS). The diagonal crosses connected by the dotted line are the Spitzer 70- jam 
derivations of Patel et al. (2013) (in the lower z-bin we report their LFs derived both at 0.0 < z < 0.2 and at 0.2 < z < 0.4). The cyan open triangles show 
the PEP 60- jam LF from the SDP data in GOODS-N by Gruppioni et al. (2010) (the redshift bins are not exactly the same). The grey dashed line shows the 
Saunders et al. (1990) local LF. 


space density of X-ray (see Brusa et al. 2009; Civano et al. 2011) 
and optically selected AGN (Richards et al. 2006), and of sub-mm 
galaxies (Wall, Pope & Scott 2008), as well as in the HerMES total 
IR LF - though less evident (from 250- jam data; Vaccari et al., in 
preparation; see Fig. 7). Moreover, our result is in agreement with 
the recent finding of Smit et al. (2012) that the characteristic value of 
the galaxy SFR exhibits a substantial, linear decrease as a function 
of redshift from z ~ 2 to z ~ 8. 

In Fig. 9, we show the values of L* and O* at different redshifts, 
with the best curve [oc (1 + z) k ] fit to the data points. The values of 
the curve slopes and of the redshifts corresponding to evolutionary 
breaks have been chosen to be those which minimize the x 2 of the 
fit with two power laws. We find a significant variation of L* with 
z, which increases as (1 + z ) 3 - 55±01 ° up to z ~ 1.85, and as (1 + 
z )i.62±o.5i | 3 e t ween z ~ 1.85 and z ~ 4. The variation of with 
redshift starts with a slow decrease as (1 + z )-°- 57±0 - 22 up to z ~ 
1.1, followed by a rapid decrease oc (1 + z) -3 ' 92 ± 0-34 at z > 1.1 and 
up to z ~ 4. 

Previous estimations of the evolution of L* and O* in the past 
(i.e. Caputi et al. 2007, Bethermin et al. 2011 and Marsden et al. 
2011) discussed a decrease in the density of far-IR sources between 
z — 1 and z = 2. In particular, Bethermin et al. (2011) and Marsden 
et al. (2011), by evolving a parametrized far-IR LF, explored the 
evolution required by the source counts in the parameter space. The 
results of these works (especially those of Bethermin et al. 2011, 
showing a decrease of at z > 1 and a flatter trend on the evolution 
of L* at z > 2) are close to ours, though with the source counts only 
it was not possible to constrain the evolution of IR sources at z > 2. 


3.5 Evolution of the different IR populations 

The PEP survey, given its size and its coverage in luminosity and 
redshift, allows us to go further in investigating the evolution of the 
total IR LF: it gives us the opportunity to study the evolution of 
the different galaxy classes that compose the global IR population. 
To investigate the different evolutionary paths of the various IR 
populations, we have computed the 1/V max LFs separately for the 
five galaxy classes defined by the SED-fitting analysis. In Fig. 10, 
we show the total IR LFs derived from the combined PEP samples 
for the different SED classes (coloured filled areas). The results of 
the fit to a modified Schechter function (see equation 3) for each 
population are overplotted on the data. In fact, by following the 
same procedure adopted for the global LF, a parametric fit to the 
LFs at different redshifts has been performed also for the single 
populations. The a and o parameters, for each population, have 
been estimated at the redshift where the corresponding LF is best 
sampled (not necessarily at the lowest z-bin as for the global LF). 
Subsequently, the values of a and o have been frozen at the values 
found in the ‘optimal’ redshift bin, leaving only L* and <3>* free 
to vary. In Fig. 11, analogously to Fig. 9, we show the values 
of L* ; and O* at different redshifts for the different populations, 
with the best least-squares fitting curves [a (1 + z) k ] overplotted. 
The best- fitting values of a, <j, L* and <3>* in the first redshift bin, 
together with the parameters describing the luminosity [k L , i , k L , 2 
and Zb,G oc (1 + z) kM to z = Zb,L, oc (1 + z) ku2 at z > Zb,iJ and 
density evolution [k Pf 1 , 2 and Zb, P - oc (1 + z) kpA out to z = Zb, P , 

oc (1 + z) kpa at z > Zb, p] are reported in Table 8. 
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Figure 5. Rest-frame 90- jam LF estimated independently from the 160- pm-selected samples in the four PEP fields with the l/V max method (with the same 
symbols as in Fig. 4). The diagonal crosses represent the local LF of Serjeant et al. (2004), the brown dots with error bars are the local LF of Sedgwick et al. 
(2011) in the AKARI Deep Field, while the pink open circles are the 100- pm LF derivation of Lapi et al. (201 1) in the H- ATLAS fields. The cyan open triangles 
show the PEP 90- pm LF computed by Gruppioni et al. (2010) from the SDP data in the GOODSN (not exactly the same redshift bins). 

Table 4. Combined rest-frame 60- pm LF. 


log(L 60 /LQ) 

0.0 < z < 0.4 

0.4 < z <0.8 

log(4VMpc 3 dex 
0.8 <z <1.2 1.2 <z <1.8 

') 

1.8 <z <2.5 

2.5 <z <3.5 

3.5 <z <4.5 

7. 5-7. 9 

-1.37 ±0.43 







7. 9-8. 3 

-1.54 ±0.32 







8. 3-8.7 

-1.82 ±0.20 







8.7-9. 1 

-1.75 ±0.09 







9. 1-9.5 

-2.22 ± 0.04 







9. 5-9. 9 

-2.36 ± 0.03 

-2.25 ± 0.27 






9.9-10.3 

-2.58 ± 0.02 

-2.70 ±0.12 






10.3-10.7 

-2.77 ± 0.02 

-2.46 ± 0.04 






10.7-11.1 

-3.28 ± 0.03 

-3.00 ± 0.02 

-2.70 ± 0.07 





11.1-11.5 

-4.13 ±0.09 

-3.41 ± 0.02 

-3.12 ±0.04 

-3.09 ±0.11 




11.5-11.9 

-4.91 ±0.22 

-4.26 ± 0.05 

-3.61 ± 0.02 

-3.47 ± 0.06 

-3.08 ± 0.22 



11.9-12.3 


-5.03 ±0.11 

-4.31 ±0.04 

-3.95 ± 0.03 

-3.63 ±0.10 

-3.22 ±0.31 


12.3-12.70 


-6.20 ± 0.43 

-5.31 ±0.12 

-4.73 ± 0.05 

-4.22 ± 0.06 

-4.09 ±0.15 


12.7-13.1 




-5.98 ±0.18 

-5.08 ± 0.06 

-4.76 ±0.13 

-4.45 ± 0.43 

13.1-13.5 





-6.88 ± 0.43 

-5.58 ± 0.22 

-5.39 ± 0.43 

13.5-13.9 






-6.93 ± 0.43 

-6.20 ± 0.43 


A clear result of our analysis is that the evolution derived for 
the global IR LF is indeed a combination of different evolution- 
ary paths: the far-IR population does not evolve all together ‘as a 
whole’, as it is often assumed in the literature, but is composed by 
different galaxy classes evolving differently and independently. As 
shown in Fig. 10, the normal spiral galaxy population dominates 
the LF at low-z, from the local Universe up to z ~ 0.5. Moving to 
higher redshifts, the number density of galaxies with spiral galaxy 


SEDs sharply decreases, while their luminosity continues to in- 
crease, at least up to z ~ 1 (see the <3>* and U parameter trends 
shown in Fig. 11). We note that what we observe between z ~ 0 
and z ~ 1 for the spiral SED galaxies is an increase of L* by a 
factor of ~5, and a decrease of <3>* by a factor of ~10. Since the 
two evolutions are not independent, the ‘global’ evolutionary ef- 
fect results from the combination of the two. By fixing the spiral 
LF at a given volume density value, we find that the logarithm of 
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Table 5. Combined rest-frame 90- jam LF. 


log(L 90 /LQ) 

0.0 < z <0.4 

0.4 < z <0.8 

log(<J>/Mpc 3 dex 1 
0.8 <z <1.2 1.2 <z <1.8 

) 

1.8 <z <2.5 

2.5 <z <3.5 

3.5 <z <4.5 

7.9-8. 3 

-1.09 ±0.31 







8. 3-8. 7 

-1.68 ±0.21 







8.7-9. 1 

-2.00 ±0.15 







9. 1-9.5 

-2.21 ± 0.07 







9.5-9. 9 

-2.31 ±0.04 







9.9-10.3 

-2.51 ±0.03 

-2.06 ± 0.24 






10.3-10.7 

-2.63 ± 0.02 

-2.57 ± 0.05 

-2.73 ±0.18 





10.7-11.1 

-3.04 ± 0.03 

-2.98 ± 0.02 

-2.87 ± 0.08 

-2.65 ± 0.22 




11.1-11.5 

-3.91 ±0.07 

-3.31 ±0.02 

-3.02 ± 0.05 

-3.06 ±0.12 

-3.34 ± 0.20 



11.5-11.9 

-4.61 ±0.15 

-4.10 ±0.04 

-3.65 ± 0.02 

-3.46 ± 0.05 

-3.36 ±0.10 

-3.34 ±0.31 


11.9-12.3 


-5.06 ±0.12 

-4.37 ± 0.04 

-4.12 ±0.03 

-3.80 ± 0.08 

-4.04 ±0.15 

-4.18 ±0.31 

12.3-12.7 



-5.76 ±0.19 

-5.02 ± 0.06 

-4.26 ±0.10 

-4.55 ±0.14 

-5.39 ± 0.43 

12.7-13.1 




-6.47 ±0.31 

-5.50 ±0.10 

-5.13 ±0.11 

-5.86 ±0.19 

13.1-13.5 






-6.61 ±0.43 
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Figure 6. Total IR LF estimated independently from the 160- pm samples in the four PEP fields (with the same symbols as in Fig. 4). The pink and sky-blue 
shaded areas represent the range of values derived with 20 iterations by allowing a change in photo-z in the GOODS-S and COSMOS fields, respectively. 


the luminosity corresponding to that value increases by a factor of 
~2.5 between z = 0 and z = 1, in good agreement with previous re- 
sults, either empirical (for morphologically classified discy galaxies; 
Scarlata et al. 2007) or theoretical (from chemical evolution mod- 
els of Milky Way-like galaxies; Colavitti, Matteucci & Murante 
2008). 


Over the whole redshift range 0.5—3, the ‘global’ LF is domi- 
nated by the SF-AGN population. The number density of SF-AGN 
is nearly constant from the local Universe up to z ~ 1-1.5, show- 
ing a slight decrease at higher redshifts, while their luminosities 
show positive evolution up to the highest redshifts (z ~ 3.5-4). 
From Fig. 10 we note a sort of bimodality in the SF-AGN LFs 
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Figure 7. Total IR LF estimated by combining the data from the four PEP fields using the Avni & Bahcall (1980) method (black filled circles). The grey 
filled area shows the uncertainty locus obtained by combining the Poissonian error with the photometric redshift uncertainty derived through Monte Carlo 
simulations. The black solid line represents our best fit to the PEP data in the different redshift bins, corresponding to the parameters reported in Table 7. Other 
results from the literature are plotted for comparison (diagonal crosses connected by a grey dashed line, LLF of Sanders et al. 2003; magenta filled stars, Le 
Floc’h et al. 2005; orange filled squares, Chapman et al. 2005; orange dot-dashed line, Caputi et al. 2007; blue open triangles, Rodighiero et al. 2010a; green 
dashed line, Magnelli et al. 2009, 2011; green open circles, Magnelli et al. 2011). Note that the Magnelli et al. (2011) data correspond to slightly different 
redshift bins: 1.3 < z < 1.8 and 1.8 < z < 2.3, respectively. We have therefore plotted the data points corresponding to the first redshift bin in our 1.2 < z < 
1.7 panel and the data point corresponding to the second redshift bin in both our 1.7 < z < 2.0 and 2.0 < z < 2.5 panels. The red open squares are the total IR 
LFs derived by Vaccari et al. (in preparation) from the HerMES 250- jam-selected COSMOS sample. 


Table 6. PEP total IR LF. 


log(Li R /LQ ) 

0.0 <z <0.3 

0.3 < z < 0.45 

log(<t>/Mpc 3 dex 1 ) 

0.45 <z <0.6 0.6 <z <0.8 0.8<z<1.0 1.0<z<1.2 1.2 <z <1.7 1.7 <z <2.0 2.0 <z <2.5 2.5 < z < 3.0 3.0 <z <4.2 

8. 5-9.0 

-2.21 ± 0.43 



9.0-9. 5 

-2.18 ±0.09 



9.5-10.0 

-2.28 ± 0.04 



10.0-10.5 

-2.50 ± 0.03 

-2.37 ±0.11 


10.5-11.0 

-2.71 ±0.02 

-2.64 ± 0.04 

-2.61 ±0.08 -2.27 ±0.15 

11.0-11.5 

-3.49 ± 0.06 

-3.12 ±0.03 

-2.99 ±0.04 -2.89 ±0.05 -3.09 ± 0.08 -2.80 ± 0.09 -2.93 ±0.18 

11.5-12.0 

-4.79 ± 0.25 

-4.29 ±0.10 

-3.89 ±0.05 -3.53 ±0.03 -3.24 ± 0.04 -3.17 ±0.06 -3.29 ± 0.06 -3.76 ±0.13 

12.0-12.5 

-5.27 ± 0.43 

-5.58 ± 0.43 

-5.50 ±0.31 -4.75 ±0.09 -4.23 ± 0.05 -4.00 ± 0.03 -3.81 ± 0.03 -3.96 ±0.11 -3.53 ± 0.08 -3.75 ± 0.21 

12.5-13.0 



-5.79 ±0.31 -5.74 ±0.25 -5.18 ±0.12 -4.85 ± 0.05 -4.42 ± 0.04 -4.40 ± 0.04 -4.15 ±0.11 -4.65 ±0.14 

13.0-13.5 



-6.48 ±0.31 -6.01 ±0.22 -5.79 ±0.13 -5.11 ±0.07 -5.75 ±0.13 

13.5-14.0 



-6.54 ±0.31 -7.18 ±0.43 
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Table 7. Parameter values describing the curve fitted to the total IR LF. 


Redshift range 

a 

o 

logio(L*/LQ) logi 0 (O*/Mpc 3 dex 1 ) 

0.0 < z < 0.3 

1.15 ±0.07 

0.52 ± 0.05 

10.12 ±0.16 

-2.29 ± 0.06 

0.3 < z < 0.45 

l.2 a 

0.5 a 

10.41 ± 0.03 

-2.31 ±0.03 

0.45 < z < 0.6 

1.2 a 

0.5 a 

10.55 ± 0.03 

-2.35 ± 0.05 

0.6 < z < 0.8 

1.2 a 

0.5 a 

10.71 ± 0.03 

-2.35 ± 0.06 

0.8 <z< 1.0 

i.2 a 

0.5 a 

10.97 ± 0.04 

-2.40 ± 0.05 

1.0 <z< 1.2 

1.2 a 

0.5 a 

11.13 ±0.04 

-2.43 ± 0.04 

1.2 <z< 1.7 

\.2 a 

0.5 a 

11.37 ±0.03 

-2.70 ± 0.04 

1.7 < z < 2.0 

1.2 a 

0.5 a 

11.50 ±0.03 

-3.00 ± 0.03 

2.0 < z < 2.5 

i.2 a 

0.5 a 

11.60 ±0.03 

-3.01 ±0.11 

2.5 < z < 3.0 

1.2 a 

0.5 a 

11.92 ±0.08 

-3.27 ±0.18 

3.0 < z < 4.2 

i.2 a 

0.5 a 

11.90 ±0.16 

-3.74 ± 0.30 


a Fixed value. 
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Figure 8. Total IR LF estimated by combining the data from the four PEP fields using the Avni & Bahcall (1980) method plotted in all the different redshift 
intervals considered in this study, from z ~ 0 to z ~ 4. The different colour filled areas represent the ± 1 a (Poissonian) uncertainty regions at different redshifts. 


(at z < 0.45, where we are able to cover a larger luminosity range). 
This bimodality is indeed to be ascribed to the crossing of two 
contributions: that of the SF-AGN(Spiral) population, responsible 
for the faint-end steepness of the LFs, and that of the SF-AGN(SB) 
population, dominating the bright-end of the SF-AGN LFs and 
declining at low L IR (not reported in the figure). 

The starburst galaxy population never dominates. The redshift 
range where we observe the highest contribution from the starburst 
galaxies is at z ~ 1-2, while in the local Universe their contribution 
is almost negligible (i.e. their parameter shows an opposite trend 
with respect to that of spiral galaxies, see Fig. 11). 

The AGN 1 and AGN2 populations show a very similar evolution- 
ary trend as a function of z, both in O* and L*. These powerful AGN 
populations dominate only the very bright end of the total IR LF, 
although their number densities and luminosities keep increasing 
from the local Universe up to the higher redshifts. At z > 2.5, the 
AGN1 and AGN2 populations become as important as the SF-AGN 
one, with the total IR LF of PACS- selected sources in the red- 
shift range 2.5-4 being totally dominated by objects containing an 
AGN. 


3.6 Total IR LF in mass and sSFR bins 

3. 6. 1 Stellar masses and SFR from SED fitting 

The wealth of multiwavelength data available in the cosmological 
fields included in our work allow us to perform a detailed SED 
fitting of all sources, in order to derive their most relevant phys- 
ical parameters (e.g. stellar masses). To derive stellar masses, we 
have fitted the broad-band spectra of our sources using a modified 
version of magphys (da Cunha, Chariot & Elbaz 2008), which is a 
code describing the SEDs using a combination of stellar light and 
emission from dust heated by stellar populations. In particular, the 
magphys software simultaneously fits the broad-band UY-to-far-IR 
observed SED of each object, ensuring an energy balance between 
the absorbed UV light and that re-emitted in the far-IR regime. The 
main assumptions are that the energy re-radiated by dust is equal 
to that absorbed, and that starlight is the only significant source of 
dust heating. We refer to Da Cunha et al. (2008) for a thorough 
formal description of how galaxy SEDs are build. At each source’s 
redshift, the code chooses among different combinations of SF his- 
tories, metallicities and dust contents, associating a wide range of 
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z 


Figure 9. Evolution of L* and <E>* as function of z (L* oc (1 + 
Z )3.55±0.10 at z < 1.85, oc (1 + z) L62±a51 a t z > 1.85; 
d>* oc (1 d- z )-°- 57±0 - 22 at z < 1.1, oc (1 + z )- 3 - 92±0 - 34 at 

z >1.1. 


optical models to a wide range of infrared spectra and comparing to 
observed photometry, seeking for x 2 minimization. Each SF history 
is parametrized in terms of an underlying continuous model with 
exponentially declining SFR, on top of which are superimposed 
random bursts (see Da Cunha et al. 2008, 2010). We note that, al- 
though the MAGPHYS assumption of exponentially declining SFR 
might not be the best to reproduce the SFR history of z > 1.5 star- 
forming galaxies (i.e. exponentially increasing or increasing SFR 
would be better choices, as widely discussed by Maraston et al. 
2010 and Reddy et al. 2012), in our specific case it does not affect 
the results. In fact, we do not use the MAGPHYS derived SFRs, 
but we compute them by integrating the best-fitting SED (resulting 
from Fe Phare). The models are distributed uniformly in metal- 
licity between 0.2 and 2 times the solar. Since the magphys code 
assumes that starlight is the only significant source of dust heating, 
thus ignoring the presence of a possible AGN component, Berta 
et al. (2013) have developed a modified version of the magphys 
code by adding a torus component to the modelled SED emission, 
combining the Da Cunha et al. (2008) original code with the Fritz, 
Franceschini & Hatziminaoglou (2006) AGN torus library (see also 
Feltre et al. 2012). The spectral fitting is performed by comparing 
the observed SED of our galaxies to every model in the generated 
library, at the corresponding redshift. A/ 2 minimization provides 



9 1011121314 9 1011121314 
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Figure 10. Total IR LF estimated with the 1 / V max method by combining the data from the four PEP fields for the different populations (their drier uncertainty 
regions are shown as coloured filled areas: green for spirals; cyan for starbursts; red for SF-AGN; magenta for AGN2 and blue for AGN1), compared to the total 
LF (drier, grey filled area, same as in Fig. 7). The best-fitting modified Schechter functions are also plotted, extrapolated to fainter and brighter luminosities 
than covered by the data (with the black curve being for the total LF and the same colours used for the filled areas as for the single populations). 


Downloaded from http://mnras.oxfordjournals.org/ at NASA Goddard Space Flight Ctr on June 18, 2013 


The PEP HerMES luminosity function - I 19 



z 


Figure 11. Evolution of L* (top) and <E>* (bottom) as a function of redshift (in the form oc (1 ± z) K ) for the different IR populations. For comparison, the L* 
and 4>* evolution of the ‘global’ total IR LF (plotted in Fig. 9) are shown as grey dashed lines. 


the quality of each fit. We must point out that the mass derivation 
for unobscured AGN (i.e. AGN1) is a problematic issue, therefore 
the masses estimated for that class of objects are the most uncertain 
ones. One source of uncertainty in the mass measurement for AGN 1 
is due to the fact that, in these objects, the UV part of the SED is 
likely dominated by the AGN rather than by the host galaxy. This 
may produce an underestimate of the mass, since, if the AGN con- 
tribution is not taken into account, the data can be fitted by a bluer, 
younger and smaller mass object. On the other hand, the mid-IR 
part of the AGN SED is dominated by dust emission from the dusty 
torus heated by the central black hole. If a proper decomposition 
into a stellar and a torus component is not performed, the use of a 
pure stellar template for estimating the mass from SED fitting tends 
to reproduce the mid-IR emission with an older, redder and larger 
mass object (mass sometimes larger by a factor of 2 than those 
derived through a decomposition procedure; Santini et al. 2012). 
These two effects lead to an increase in the uncertainty of the mass 
derivation, although they might somewhat compensate their effects 
for a large sample of objects. Bongiorno et al. (2012) found that 
for the majority of COSMOS AGN1, the average ratio between the 
mass derived with and without taking into account the AGN compo- 
nent is ~1, although with a large dispersion, while for ~50 per cent 
the mass is underestimated by factors as high as 3-10 if the AGN 
contribution is not taken into account. For this reason, we obtained 
measurements of the stellar masses of our objects containing an 
AGN by means of the specific decomposition technique developed 
by Berta et al. (2013), to separate stellar and nuclear emission com- 


ponents. Examples of the results of this decomposition applied to 
SF-AGN, AGN2 and AGN1 are shown in Fig. 12. Masses of AGN 
estimated with the original magphys and with the Berta et al. (2013) 
code have been compared, showing very good agreement and small 
dispersion around the 1-1 relation. Similarly, as further check, we 
have also computed stellar masses with different code (hyperz; Bol- 
zonella, Miralles & Pello 2000) and stellar library (BC03, Bruzual 
& Chariot 2003, instead of the CB07 used as default by magphys), 
finding good agreement and no systematics, too. 

We have derived the SFRs from the total IR luminosities (es- 
timated from the SED fitting described in Section 2.1) with the 
standard Kennicutt (1998) relation (converted to Chabrier IMF), 
after subtracting the AGN contribution to L IR . We note that the total 
IR luminosity in PEP sources is usually dominated by SF, even in 
objects for which an AGN dominates the optical/near-IR/mid-IR 
part of the spectrum. 

3.6.2 LFs in different mass bins 

We compute the total IR LF for galaxies of different stellar masses: 
8.5 < log(M/M 0 ) < 10, 10 < log(M/M 0 ) < 11 and 11 < 
log(M/MQ) < 12, by means of the standard 1/V ma x formalism, 
and we show the results in Figs 13 (total IR LF in different z- 
bins) and 14 (ratio between the LF in the mass intervals and the 
total IR LF). We have compared our results with the SFR function 
[SFR converted to IR luminosity using the Kennicutt (1998) rela- 
tion] of massive (log(M/MQ >10) galaxies derived by Fontanot 


Table 8. Parameter values describing the curve fitted to the total IR LF of the different SED populations. 



a 

or log io (±*/ 

L o) 

(0.0 < Z < 0.3) 

logio(0*/ 
Mpc -3 dex -1 ) 

kh, l 

&L,2 

Zb,L 

k p , i 

k p , 2 

£b, p 

Spiral 

1.00 ±0.05 

0.50 ± 0.01 

9.78 ± 0.04 

-2.12 ±0.01 

4.49 ±0.15 

0.00 ± 0.46 

1.1 

-0.54 ±0.12 

-7.13 ±0.24 

0.53 

Starburst 

1.00 ±0.20 

0.35 ±0.10 

11.17 ±0.16 

-4.46 ± 0.06 

1.96 ±0.13 



3.79 ± 0.21 

-1.06 ±0.05 

1.1 

SF-AGN 

1.20 ±0.02 

0.40 ±0.10 

10.80 ± 0.02 

-3.20 ± 0.01 

3.17 ±0.04 



0.67 ± 0.05 

-3.17 ±0.15 

1.1 

AGN2 

1.20 ±0.20 

0.70 ± 0.20 

10.80 ± 0.20 

-5.14 ±0.17 

1.41 ±0.33 



2.65 ± 0.32 



AGN1 

1.40 ±0.30 

0.70 ± 0.20 

10.50 ± 0.20 

-5.21 ±0.11 

1.31 ±0.02 



3.00 ± 0.25 
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z= 1.22 class = SF— AGN 


z= 0.38 class=AGN2 


z= 1.25 class=AGN 1 





Figure 12. Examples of PEP AGN SEDs fitted by the magphys code (Da Cunha et al. 2008) modified by Berta et al. (2013) to add an AGN component from 
the Fritz et al. (2006) library. From left to right the result of the fit to an SF-AGN, an AGN2 and an AGN1 is shown. The stellar component unattenuated by 
dust is shown in green, while the dust-attenuated spectrum with dust IR re-emission is shown in blue. The AGN contribution is shown in red, while the black 
curve represents the total fitted spectrum, obtained from the sum of the blue and red components. Our data are represented by the black dots. 
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Figure 13. Contribution to the total IR LF from galaxies in three different mass intervals: 8.5 < log(M/MQ) <10 (green); 10 < log(M/MQ) <11 (orange) 
and 11 < log(M/MQ) < 12 (red). For comparison, the SFR function of log(M/MQ)>10 sources in the GOODS-S by Fontanot et al. 2012 (converted to total 
IR LF) is shown as black open squares. 


et al. (2012) from the GOODS-MUSIC sample at redshift 0.4 < 
z < 1.8. In the common redshift and luminosity range, we find 
an excellent agreement with our total IR LF, which is dominated 
by sources with 10 < log(M/Mo) < 11. At lower luminosities, 
not sampled by our data, the Fontanot et al. (2012) LFs are char- 
acterized by a double-peaked structure, interpreted in terms of 


the well-known bimodality in the colour(SFR)-mass diagram. As 
expected from the SFR-mass relation, the knee (L*) of our IR 
LFs in different mass bins moves to higher luminosities with in- 
creasing masses [i.e. at 0.0 < z < 0.3, log(L*/LQ) changes from 
~9.5 for log(M/Mo) = 8.5-10 sources, to ~11.3 for sources with 
log(M/M 0 ) > 11]. 
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Figure 14. Ratio between the IR LF of galaxies in three different mass 
intervals: 8.5 < log(M/MQ) < 10 (top); 10 < log(M/MQ) <11 (middle) 
and 11 < log(Af/MQ) <12 (bottom), as shown in Fig. 13, and the total 
IR LF, plotted in all the different redshift intervals considered in this study, 
from z ~ 0 to z ~ 4. The different colours represent the different redshifts, 
as explained in the plot. 


The slope of the total IR LF in each mass bin is always similar to 
(or flatter than) the ‘global’ LF (total, including all the masses). The 
lower mass galaxies dominate at lower luminosities [log(L/LQ] < 
9), while the most massive galaxies [log(Af/Mo) >11] contribute 
only at higher luminosities, even if they never dominate the LF. 
At all masses, the LF evolves with redshift, following the evolu- 
tion of the ‘global’ LF. Fig. 14 shows that the main contribution 
(>50 per cent) to the total IR LF is due to intermediate-mass ob- 
jects [log(M/MQ] = 10-11) at all redshifts and luminosities, with 
their fraction remaining almost the same from z = 0 to z ~ 4, simply 
shifting to higher luminosities. Lower mass objects [log(M/MQ) = 
8.5-10] contribute significantly only at log(L/LQ) <10, with their 
fraction just shifting to higher luminosities with redshift, but al- 
ways being below 30 per cent at z > 0.45 and log(L/LQ) >11. 
The contribution of the most massive objects [log(Af/MQ) =11- 
12] increases with IR luminosity and redshift, becoming significant 
(>50 per cent) only at z > 1.7 and log(L/LQ) > 12.5. 


3.7 Specific-SFR and the MS of star-forming galaxies 

Having computed stellar masses and SFRs for each source, we 
can check how the PACS-selected sources populate the SFR-stellar 
mass plane and the so called main sequence (MS) of star-forming 
galaxies, as a function of redshift. This relation (i.e. SFR versus 
stellar mass) has been shown to be quite tight in the local Universe 
(Peng et al. 2010, 201 1) and well established at redshift z ~ 1 (Elbaz 
et al. 2007) and up to z ~ 2 (Daddi et al. 2007; Rodighiero et al. 
2011) and z ~ 3 (Magdis et al. 2010), with normalization scaling 
as ~(1 + z) 2 ' 8 out to z ~ 2, as shown by Sargent et al. (2012) 
(see also Elbaz et al. 2007; Pannella et al. 2009; Rodighiero et al. 
2010b; Karim et al. 2011). At z ~ 2 and ~1.5, we assume a slope 
of 0.79 for the MS in the SFR versus stellar-mass plane (according 
to Rodighiero et al. 2011 and Sargent et al. 2012), while at z ~ 1 


we assume a slope of 0.9, as found by Elbaz et al. (2007), and we 
limit our investigation to the redshift range 0.8 < z < 2.2. 

By combining UV and far-IR data, Rodighiero et al. (2011) re- 
evaluated the locus of the MS at z ~ 2, showing that objects lying a 
factor of 4 above the MS (in SFR) can be considered as outliers with 
respect to the average locus where smoothly star-forming galaxies 
spend most of their lives in a secular and steady regime. In that 
work, off- sequence sources (characterized by very high sSFRs) are 
assumed to be in a starburst mode, and are found to contribute only 
2 per cent of mass-selected star-forming galaxies and to account for 
only 10 per cent of the cosmic SFR density at z ~ 2 (Rodighiero 
etal. 2011). 

In order to check what kind of objects we could classify as on- 
MS and off-MS sources in our IR sample compared to previous 
findings, either based or not based on IR surveys (e.g. Rodighiero 
et al. 2011; Sargent et al. 2012), we have split our sample into 
off-MS and on-MS. For consistency with previous studies, we have 
applied the same criterion as Rodighiero et al. (201 1) (0.6 dex above 
the MS) over the whole 0.8 < z < 2.2 redshift range, by using 
as a reference MS the relation found by Rodighiero et al. (2011) 
at z ~ 2, scaled as described above at z ~ 1 .5, and the relation found 
by Elbaz et al. (2007) at z ~ 1. In Fig. 15, we show the SFR versus 
stellar mass distributions in three redshift bins (0.8 < z < 1.25, 
1.25 < z < 1.8 and 1.8 < z < 2.2), for the PACS sources included 
in the computation of the luminosity functions presented in this 
work. The colour code marks the different SED-classes to which 
each source belongs. We also report the typical loci of the MS at the 
various redshifts [scaling as (1+z) 2 ' 8 , as mentioned above]. Details 
are given in the caption to the figure. The typical far-IR selection bias 
(PACS -Herschel in this case) appears as an approximate horizontal 
SFR cut (Rodighiero et al. 2011, Wuyts et al. 2011), shown as thin 
dotted line in Fig. 15. We note that the trends of mid/far- IR SEDs 
with offset from the main sequence observed in Fig. 15 (and widely 
discussed in Section 5) are in good agreement with the results of 
Elbaz et al. (2011) and Nordon et al. (2012). 

With the selection presented in Fig. 15 (sources qualify as ‘off- 
MS’ if they lie more than 0.6 dex above the observed SFR-stellar 
mass relation), we can compute the contribution of off-MS (also 
called ‘starburst’ in the literature) and on-MS (‘steady star-formers’ 
in the literature) sources to the total IR LFs. This is presented in 
Fig. 16, where the total IR LFs of on-MS and off-MS sources have 
been computed independently for the three redshift bins. The pink 
and yellow filled areas correspond to the ±lcr uncertainty regions 
of the total IR LFs estimated for the on-MS and off-MS populations, 
respectively, while the black filled area marks the global population. 
Our results are compared with the recent predictions by Sargent 
et al. (2012), a non-parametric approach that is based on three basic 
observables: the redshift evolution of the stellar MF for star-forming 
galaxies; the evolution of the sSFR of MS galaxies and a double- 
Gaussian decomposition of the sSFR distribution at fixed stellar 
mass into a contribution (assumed redshift- and mass-invariant) 
from MS and off-MS (i.e. starburst) activity. The evolution of the 
two populations found both in our data and the Sargent et al. (2012) 
model are very similar. Both data and predictions indicate that the 
bright-end of the total IR LF is dominated by off-MS sources. 
However, although consistent within the uncertainties, the relative 
contribution of off-MS sources seems to be stronger in the Sargent 
et al. (2012) model than observed in the present computation, where 
we find that the bright-ends of the PEP IR LFs are more similarly 
populated by MS and off-MS sources (especially at z ~ 2). This 
difference can be at least partly ascribed to the sharp cut we apply 
to separate MS from off-MS sources, while Sargent et al. (2012) 
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Figure 15. SFR versus stellar mass for our PEP 160- pm sources (green, 
spiral; cyan, starburst; red, SF-AGN; magenta, AGN2; blue, AGN1), in three 
redshift bins (from left to right): 0.8 < z < 1.25; 1.25 < z < 1.8 and 1.8 < 
z < 2.2. The relation known as MS is plotted as a solid line (from Elbaz 
et al. 2007 in the lower redshift bin, rescaled as (1 + z) 2-8 in the central 
bin and from Rodighiero et al. 2011 at z ~ 2), while the dashed line shows 
the same relation 0.6 dex higher, indicating the separation between MS and 
above MS objects adopted by Rodighiero et al. (2011). The horizontal dotted 
lines show the nominal SFR limits of the GOODS-S (lower) and COSMOS 
(upper) samples in the different redshift intervals. The orange open squares 
and the dark-green open diamonds mark the two sub-classes of SF-AGN 
galaxies: the SF-AGN(SB) and SF-AGN(Spiral), respectively. 


model the off-MS and on-MS populations with two continuous log- 
normal distributions centred at 0.6 dex above the MS and exactly 
on the MS respectively, of which the one describing the ‘starburst’ 
population has wings that extend into our on-MS selection region 
(hence attributing more sources to the starburst category than are 
selected in our off-MS class). Better agreement between our data 
and the predictions of Sargent et al. (2012) is found at the faint- 
end of the LFs that appear to be completely dominated by the 
‘normal’ MS galaxies at all redshifts (although the total and relative 
contributions at log(Li R /LQ) < 12 are slightly lower in the data 
than predicted by Sargent et al. 2012). Good agreement is also 
found with respect to the evolution of the cross-over luminosity (i.e. 


where the contributions from on-MS and off-MS sources are equal); 
in the model of Sargent et al. (2012), the cross-over luminosity 
simply shifts to higher luminosities and lower densities (according 
to the luminosity and density evolution considered). Note that the 
model assumptions rely on results from different surveys, selected at 
different wavelengths, complete in mass and with good sampling of 
the MS. On the other hand, our selection is in SFR and our sources do 
not follow any clear sequence in stellar mass-SFR. These different 
selection effects are highly likely to lead to also some differences 
between our LFs and the predictions of Sargent et al. (2012). 

To quantify the relative contribution of the two populations, our 
observed data have been fitted with a modified Schechter function, 
in order to integrate them and compute their comoving number and 
luminosity densities as functions of redshift (see the next section). 


4 NUMBER DENSITY AND IR LUMINOSITY 
DENSITY 

We derive the evolution of the comoving number and luminosity 
density (total, see Fig 18) of the PEP sources, either belonging to 
the different SED classes (Fig 18, left), to the on-MS and off-MS 
categories (middle) and to the different mass intervals (right), by 
integrating the total IR LF in the different redshift bins from z ~ 0 
to z ~ 4. To compute the number (and IR luminosity) density, we 
integrate the Schechter functions that best reproduce the different 
populations/mass/sSFR-classes, down to log(L/LQ) = 8. We note 
that here we consider lower limits the number and luminosity den- 
sities at 3.0 < z < 4.2, since our LF estimate in that redshift bin is 
likely to be incomplete, as discussed in Section 3.3. We find that the 
number density of the whole IR population is nearly constant in the 
z = 0-1.2 redshift range (slightly increasing from z ~ 0 to z ~ 0.5), 
decreasing at z > 1.2 (see top panels of Fig. 18). When decompos- 
ing the number density according to the different SED classes, we 
observe that normal spiral galaxies dominate the local density, with 
a smaller contribution also from the SF-AGN population and a neg- 
ligible one from starburst, AGN1 and AGN2. The space density of 
spiral galaxies decreases rapidly at z > 0.5, while that of SF-AGN 
keeps nearly constant at 0.5 < z < 2.5, largely dominating in that 
redshift range. Starburst galaxies never dominate, while the number 
density of the bright AGN (both AGN1 and AGN2) increases with 
redshift, from ~10 -4 Mpc -3 at z ~ 0 to ~l-2 x 10 -3 Mpc~ 3 at z 
~ 3. At higher redshifts the AGN population largely dominates the 
number density. 

If the overall contribution to the IR luminosity density (/?ir) from 
the AGN components of galaxies is small, pi R can be considered 
as a proxy of the SFR density (psfr)- As a further check, we have 
therefore studied the evolution of the SF-AGN population (which 
dominates the distribution of sources) by dividing this class into 
SF-AGN(SB) and SF-AGN(Spiral) sub-classes and studying their 
evolution separately. Indeed, we have found different evolutionary 
paths for the two populations, the former dominating at higher red- 
shifts and showing a behaviour similar to that of AGN-dominated 
sources (e.g. AGN1 and AGN2), the latter dominating at interme- 
diate redshifts (between z ~ 1 and 2), rising sharply from z ~ 2 
towards the lower redshifts and decreasing as the spiral population 
rises at z < 1. These evolutionary trends, in terms of number and 
luminosity density, have been reported in Fig. 18 as orange solid 
[SF-AGN(SB)] and dark-green dashed [SF-AGN(Spiral)] curves. 

Galaxies following the SFR-mass relation are always dominant 
over the off-MS population, at all redshifts (although their space 
density decreases with increasing z, as well as the ‘global’ number 
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Figure 16. Contribution of MS (±lcx uncertainty region: pink shaded area) and off-MS (±1 o uncertainty region: yellow shaded area) galaxies to the total IR 
LF (black filled dots and area) in three different redshift bins. For comparison, the recent predictions of Sargent et al. (2012) in similar z-bins are shown as 
grey (total, in background), light grey (MS) and dark grey (off-MS) filled regions. 


density), while the number density of the latter population remains 
nearly constant between z ~ 0.8 and z ~ 2.2. 

In all the mass bins, the trends with redshift of the galaxy num- 
ber densities are similar to the ‘global’ one, decreasing at higher 
redshifts, although with slightly different slopes for the different 
mass intervals. The number densities of low-mass galaxies [8.5 < 
log(M/Mo) < 10], reported in the top-right panel of Fig. 18, have 
been computed by integrating the best-fitting modified Schechter 
function only to z ~ 2, since data were not enough to derive reliable 
fits at higher redshifts. To this redshift, these sources outnumber 
the higher mass ones, although they fall steeply above z ~ 1, when 
they reach about the same volume density of higher mass galaxies 
[10 < log(M/MQ) < 11]. Massive objects [log(M/MQ) >11] never 
dominate (always below 5 per cent) the total number density. 

The total IR LF allows a direct estimate of the total comoving IR 
luminosity density (pi R ) as a function of z, which is a crucial tool 
for understanding galaxy formation and evolution. Although pi R 
can be converted to an SFR density (psfr) under the assumption 
that the SFR and Li R quantities are connected by the Kennicutt 
(1998) relation, before doing that we must be sure that the total 
IR luminosity is produced uniquely by SF, without contamination 
from an AGN. The SED decomposition and separation into AGN 
and SF contributions show a negligible contribution to Li R (< 10 per 
cent) from the AGN in most of the SF-AGN, and an SF component 
dominating the far-IR even in the majority of more powerful AGN 
(AGN1 and AGN2). Here, we prefer to speak in terms of pi R rather 
than of Psfr, since, especially at high redshift - where the AGN- 
dominated sources are more numerous - the conversion of pi R could 
represent only an upper limit to Psfr- Note, however, that since this 
population is never dominant in our IR survey, we do not expect 
that contamination related to accretion activity occurring in these 
objects (mainly at high-z) can significantly affect the results in terms 
of Psfr- 

In Fig. 17, we show pi R estimated from our total IR LF and com- 
pare it with results obtained from previous IR surveys (Le Floc’h 
et al. 2005; Caputi et al. 2007; Rodighiero et al. 2010a; Magnelli 
et al. 2011). In the common redshift intervals (0 < z < 2-2.5), we 
find very close agreement with previous results based on IR data, 
especially with the Magnelli et al. (201 1) derivation. As well as pre- 
vious findings, pi R from PEP shows the rapid rise from z ~ 0 to z ~ 
1, followed by a flattening at higher redshifts. The indications from 


i i i i 
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Figure 17. Redshift evolution of the total IR luminosity density (pi R , ob- 
tained by integrating the modified Schechter functions that best reproduce 
the total IR LF down to log(L/LQ) = 8) to z = 4. The results of integrating 
the best-fitting curve for our observed total IR LF in each z-bin are shown as 
black filled circles (the grey filled area represents the ± io uncertainty locus) 
and compared with estimates from previous mid-IR surveys (magenta filled 
area, Le Floc’h et al. 2005; orange filled triangles, Caputi et al. 2007; blue 
open triangles, Rodighiero et al. 2010a and green open circles, Magnelli 
et al. 2011). The upward pointing arrow in the highest-z bin means that, 
due to the large fraction of photometric redshifts and the fact that the PEP 
selection might miss high-z sources, our 3.0 < z < 4.2 pi R estimate is likely 
to be a lower limit. 

our survey are that the intermediate-redshift flattening is followed 
by a high-redshift decline, which starts around z ~ 3. From our 
data, pi R evolves as (1 + z ) 3 0±0 - 2 up to z ~ 1.1, as (1 + z )-°- 3±01 
from z ~ 1.1 to z ~ 2.8, and then as (1 + z )- 6 0±0 - 9 up to z ~ 4. 

In the bottom panels of Fig. 18, we plot the different contribu- 
tions to pi R from the different SED populations (left), from the 
on-MS and off-MS sources (middle) and from the different mass 
intervals. We notice a predominance of spiral- SED galaxies only 
at low redshifts (z < 0.5-0.6), when SF-AGN begin to dominate 
Pir up to z ~ 2.5. The starburst SED galaxies are never the preva- 
lent population, although their contribution to pi R increases rapidly 
from the local Universe to z ~ 1, then keeps nearly constant to 
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Figure 18. Top: evolution of the comoving number density of PEP sources up to z ~ 4 (black filled circles with error bars within the ±lcr uncertainty region, 
represented by the grey filled area). The upward pointing arrow in the highest- z bin means that our 3.0 < z < 4.2 estimates are likely to be lower limits. Bottom: 
redshift evolution of the total IR luminosity density to z = 4. To compute the number and IR luminosity density, we integrate the modified Schechter functions 
that best reproduce the different populations/mass/sSFR-classes, down to log(L/LQ) = 8. The black filled circles and the grey dashed area in all the three 
panels represent our PEP derived pir and its ±lcr uncertainty region, as shown in Fig. 17. In the left-hand panels, we show the number (top) and luminosity 
density (bottom) of the different IR populations (green filled area, spiral; cyan, starburst; red, SF-AGN; magenta, AGN2; and blue, AGN1). The contribution of 
SF-AGN sources, sub-divided on the basis of their SED resemblance to spiral or starburst templates, is shown by the dark-green dashed [SF-AGN(Spiral)] and 
orange solid [SF-AGN(SB)] lines, respectively. In the middle panels, we show the uncertainty regions of the relative contribution to number and luminosity 
density of the sources on- and off- the SFR-stellar mass MS (Daddi et al. 2007; Elbaz et al. 2007), as pink and yellow filled areas, respectively. Our derivations 
have been compared to those of Bethermin et al. (2012) for on-MS and off-MS sources (converted from psfr to pir using the Kennicutt 1998 relation), which 
are represented by the purple and orange dashed lines, respectively. In the right-hand panels, we show the relative contribution to the number and luminosity 
density of sources with different masses (green, 8.5 < log(A//MQ) < 10; orange, 10 < log(A//MQ) <11; and red, 11 < log(M/MQ) < 12). For comparison, 
in the bottom-right panel, we plot also the results of Santini et al. (2009) in the GOODS-S field in similar mass intervals (light-orange triangles and dashed 
line: 9.77 < log(A/ZMQ) < 10.77; pink triangles and dashed line: log(A//MQ) > 10.77). 


z ~ 2.5, to decrease at higher redshifts. The SF-AGN(SB) and SF- 
AGN(Spiral) contributions to pi R show opposite trends, with the 
former sharply increasing towards the higher redshifts (dominating 
at z > 2), and the latter prevailing between z ~ 1 and ~2, then 
dropping at higher redshifts. AGN 1 and AGN2 start dominating the 
IR luminosity density at z > 2.5, with their pi R always rising from 
z ~ 0, and then remaining almost constant (or slightly decreasing) 
towards the higher redshifts. Note how these two populations of 
AGN evolve similarly, as an indication of their same intrinsic na- 
ture and of the absence of any significant bias in the far-IR selection. 

The contributions to pi R of MS and off-MS sources (pink and 
yellow filled regions, respectively, in the bottom-middle panel of 
Fig. 18) keep nearly constant between z ~ 0.8 and z ~ 2.2. In par- 
ticular, the off-MS sources contribution keeps around 20 per cent of 
the total Pir over the whole 0.8 < z < 2.2 redshift interval, showing 
no significant signs of increase (or decrease). Our results strengthen 


the role of MS sources (pink shaded regions) in the build up of the 
stellar mass in galaxies, at all cosmic epochs, with evidence that 
their role is even increasing from z ~ 2 to z ~ 1: their number 
density changes by a factor of ~2, while their luminosity /SFR den- 
sity remains nearly constant. The importance of the off-MS sources 
does not show any significant signs of decreasing at lower z, their 
relative (with respect to the total) number (luminosity) density pass- 
ing from ~9 per cent (22 per cent) at z ~ 2 to ~6 per cent (19 per 
cent) at z ~ 1. These fractions are relatively different from those 
found by Rodighiero et al. (2011) (off-MS galaxies represent only 
2 per cent of mass-selected star-forming galaxies and account for 
only 10 per cent of the cosmic SFR density at z ~ 2). However, 
we must note that Rodighiero et al. (2011) for their analysis com- 
bined far-IR-selected (i.e. SFR-selected) and near-IR- selected (i.e. 
M*-selected) star-forming samples, well defining the MS, while 
with our data (SFR-selected only), we are not able to observe any 
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correlation between SFR and stellar mass (see Fig. 15) and barely 
detect MS objects at z ~ 2. 

Our results have been compared to the SFR densities (converted 
to Pir using the Kennicutt 1998 relation) for on-MS and off-MS 
sources based on the Sargent et al. (2012) model recently derived 
by Bethermin et al. (2012). These are shown in Fig. 18 as purple 
and orange dashed lines, respectively. The predicted off-MS psfr 
agrees well with our estimate in the common redshift range, while 
the predicted one for MS sources is higher than that derived from 
our data (especially at z ~ 1.5). As already discussed regarding the 
comparison with the Sargent et al. (2012) model, this discrepancy is 
likely to be ascribed to our difficulty to extrapolate the MS to low- 
SFR values and to the different selection criteria used to separate 
MS from off-MS sources. 

In the bottom-right panels of Fig. 18, we show the contribution of 
the different mass populations to the luminosity density as a function 
of redshift. Although we detect a similar steep increase of pi R versus 
redshift at z < 1 in both low- and high-mass galaxies, the evolution 
in Pir of galaxies with different masses is very different, reflecting 
the downsizing scenario, with pi R peaking at higher redshift with 
increasing mass. Indeed, the IR luminosity density of intermediate- 
mass objects [log(M/MQ) = 10-11] always dominates, increasing 
up to z ~ 1, then remaining nearly constant at higher redshifts (at 
least up to z ~ 2.8). The IR luminosity density of most massive 
objects increases even more rapidly with redshift (at z = 2 it was 
approximately five times higher than today) and continues to grow 
up to z = 3, where their contribution to pi R is ~30 per cent of the total 
and close to that of intermediate-mass objects (which contribute 
~60 per cent at z > 3). We compare our results with those of 
Santini et al. (2009), plotted in the figure as thick dashed lines [light 
orange: 9.77 < log(M/MQ) < 10.77; pink: log(M/MQ) > 10.77; 
with a Chabrier IMF used to determine our masses], showing very 
similar trends and values for both intermediate- and higher mass 
galaxies. Our analysis of high-mass galaxies extends up to z ~ 4, 
finding that for the most massive galaxies pi R continues to rise even 
at z > 2, with an apparent peak at z = 3. This result confirms that 
the formation epoch of galaxies proceeded from high- to low-mass 
systems. 

The values of pi R in the different redshift intervals, either the 
total ones or the contributions from the different classes (SED, 
mass, sSFR), are reported in Table 9. 

5 DISCUSSION 

In the previous sections, we have discussed the different evolu- 
tionary behaviour of different classes of sources, either divided by 
SED-type, mass or sSFR. In this section, we will try to understand 
‘who is who’, discussing which populations are mainly on-MS or 
off-MS, which have the larger (smaller) masses, and how and if the 
relative contributions of these populations vary with redshift. From 
Fig. 15, we note that the off-MS sources are dominated by galaxies 
with AGN-type SEDs. In the lower redshift bin considered, MS 
sources are mainly spirals and SF-AGN(Spiral), with a small tail 
of these populations contributing also to the off-MS. AGN1 and 
AGN2 are prevalently 0.6 dex above the MS (off-MS). Although 
the mass estimate of AGN1 suffers from the largest uncertainties, 
our result suggests that the off-MS population is largely constituted 
by sources containing an AGN, with a higher fraction of AGN- 
dominated objects at higher z. This should indicate that the major 
merger episode likely associated with what is considered a different 
mode of SF (e.g. Wuyts et al. 2011 finds that off-MS galaxies are 
likely mergers from their Sersic index) could also trigger an intense 




B 

o 

Ph 

w 

Oh 


04 

00 



ON 

O' 

04 

in 




o 




4 




or 

’ — 1 

OO 

o- 



o- 

’ — 1 




V 

4 



d 

d 

d 

d 



d 

d 




IV 

-H 



4 

4 

4 

4 



4 

4 




VI 

o 



ON 

ON 

00 

04 



00 

or 





in 



04 

CO 

04 

ON 



q 

CO 




co 

04 




d 

’ —t 

d 




d 




q 

CO 


CO 

00 


ON 

04 

04 


ON 

or 




CO 

CO 



in 

CO 

in 

NO 

CO 


q 

04 




V 

CO 


d 

d 

d 

d 

d 

d 


04 

04 




tv 

4 


4 

4 

4 

4 

4 

4 


4 

4 




VI 

o 


o 



o- 

in 

CO 



04 




in 

04 

00 


OT 

q 

04 

o 

CO 

O' 


or 

00 




4 


d 

CO 


04 

04 

d 


4 

04 




in 

in 


04 

o 

04 

OO 

of 

in 



NO 




04* 



CO 

NO 

04 

1 — 1 

O' 

in 


NO 

or 




V 

4 


d 

d 

d 

d 

d 

d 


d 

d 




tv 

41 


4 

4 

4 

4 

4 

4 


4 

4 




VI 

^ o 


in 



in 

CO 

oo 


oo 

CO 





o 


in 


or 

04 

00 

00 



04 




04 

4 


1 

in 

d 

1 

CO 

d 


CO 

,—l 




q 



04 

CO 


oo 

04 

oo 



04 




c4 

or 



04 


o 


o 


or 

04 




V 

d 


d 

d 

d 

d 

d 

d 


d 

d 




IV 

4 


4 

4 

4 

4 

4 

4 


4 

4 




VI 

in 


00 

in 

o 

ON 

04 

00 


00 

04 





00 



i> 

NO 

in 


CO 


04 

OO 




- 

4 


d 

04 

d 

d 


d 


04 

d 





CO 

CO 

o~ 

ON 

o 

NO 

NO 

ON 

04 

in 

o 




' — 1 


o 

o 

o 


o 

o 

04 

in 

oo 

04 




V 

d 

d 

d 

d 

d 

d 

d 

d 

d 

4 

d 




IV 

41 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 




VI 

ON 

04 

04 

oo 

in 

or 


ON 

04 

oo 

in 




04 

CO 

CO 


04 

04 

or 

00 

q 

q 

1 — 1 

oo 




00 

d 

d 

04 

d 

d 

d 

04 


4 

d 




04 

in 

o 

CO 

NO 

04 

o 

in 

00 


o 





,—l 

on 


1 — 1 

CO 

o 



CO 

ON 

ON 

CO 




V 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 




cC ** 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 




1 VI 

04 

or 

o 


oo 

04 

in 

CO 


in 





o 

Oh O 

s - 

O 

on 

i> 

00 

o 


q 

04 

04 

q 

o 




On 

d 

d 

CO 

d 

d 


04 


4 

d 




G) 















4 p 

CO 

ON 

OT 

04 


04 

in 

NO 


04 





00 ^ 

00 

o 

o 


O 

o 

o 


O' 

in 

o 




2 v 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 




q tv 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 




& VI 

, i 

oo 

OT 

04 

CO 

in 

NO 

NO 

i 

or 

CO 




^ oo 

v o 

ON 

in 

00 

o 

o 


04 

q 

O) 

NO 




d 

no 

d 

d 

04 

d 

d 

d 


** 

04 

d 




oo 

00 

in 

CO 

04 

OT 

04 

CO 

00 

in 

in 

00 




d 

in 


o 


o 

O 

o 

1 — 1 

CO 

o- 

1 — 1 




V 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 




IV 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 




VI 

t> 

04 

oo 

00 

00 

O' 

CO 

NO 

i 

or 

i 




NO 

d 


q 

1 — 1 

04 

o 

o 

04 

in 

i — i 

in 

CO 




4 

04 

d 

04 

d 

d 

d 



CO 

d 




q 

o 

04 

04 

04 

NO 





'xf 

o 

04 

04 

CO 

in 


CO 


O 

O 


o 

o 

o 


or 


04 

q 

NO 

V 

d 

d 

d 

d 


d 

d 

d 

d 

d 

d 

V 

in 

d 

IV 

4 

4 

4 

4 


4 

4 

4 

4 

4 

4 

IV 

4 

4 

VI 

(N 

00 

o 

04 


oT 

o 

O' 


or 

CO 

v 

04 

i 

in 

00 

q 


in 


O 


oo 

ON 


04 

oo 


q 

or 

d 

04 

7—1 

d 



d 

d 

d 

d 

04 

d 


4 

04 

in 

o^j- 





or 

o- 


in 

o 

NO 

CO 


NO 

o 

00 



d 

2 

o 

o 

o 

o 

o 

o 

04 


04 

So 

v 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

V 

4 

d 

IV 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

IV 

4 

4 

VI 

in 

ON 

NO 

CO 

04 

ON 

NO 

o 

or 

in 

NO 

V 

04 

ON 

04 

CO 

o 

o 


04 

o 

NO 

O' 

NO 

04 

in 

q 


CO 

oi 

4 

d 

4 

o 

O 

d 

d 

d 

4 

d 

04 

4 

q 

d 





d 

d 






1 



CO 

ON 

(N 





ON 


04 

O 




in 

04 



d 

o 

o 

o 

o 

o 

s 

o 

o 

o 

q 

2 

V 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 


4 

d 

IV 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

V 

tv 

4 

4 

VI 

no 

o 

in 

o 

04 

in 

NO 

o 


ON 

00 


i 



CO 

q 

o 

or 

O 

04 

o 

"xl- 

CO 

ON 

o 

V 

i> 

or 

P 

4 


d 

d 

d 

o 

d 

d 

d 

d 

d 

00 

in 

4 

o 






d 



o 


04 

d 












T 

in 

00 

II 

3 

T 












S-H 

ii 

ii 











pq 

‘Bh 














cn 


o o o 





4 

gs 

l 

co 

q 

Z 

% 

04 

Z 

Z 

Z 

5? 

z 

% 

s s s 




o 

’Bh 

£ 

Ph 

Q 

o 

Ph 

Ph 

w> 

(oO 

OX) 


1=1 

d 


H 

cn 


00 

< 

< 

cn 

00 

o 

o 

o 


O 

o 


Downloaded from http://mnras.oxfordjournals.org/ at NASA Goddard Space Flight Ctr on June 18, 2013 


26 


C. Gruppioni et al. 



8 9 10 11 12 



9 1011 9 1011 9 1011 

log(M/M 0 ) log(M/M 0 ) log(M/M 0 ) 


Figure 19. Mass distribution of the PEP 160- pm sources in all the four PEP fields, with different colours corresponding to different SED-classes (green, spiral; 
cyan, starburst; red, SF-AGN; magenta, AGN2; and blue, AGN1). The top panel shows the mass distribution of the different IR populations at all redshifts, 
while the bottom panels report the mass distribution in three representative redshift bins (0 < z < 0.3, 0.8 < z < 1.0 and 2.0 < z < 2.5). 


AGN activity, whose presence strongly influences the SFR within 
the host galaxy. While AGN show preferentially high sSFRs, in 
our sample they do not seem to prefer higher mass systems (see 
Fig. 19). Keeping in mind that the masses in our AGN1 objects are 
affected by large uncertainties, we can interpret this as due to the 
fact that ours is a selection in SFR. Therefore, our AGN are still in 
an intense SF phase, where they are still accreting their stellar mass 
by actively forming stars. Due to the far-IR selection, in fact, we 
miss the further phase of quiescent, more massive population, with 
the SF totally quenched and all the stellar mass already formed and 
in place. Moreover, our result is not in conflict with the consensus in 
the literature (i.e. that X-ray AGN above a certain X-ray luminosity 
do prefer massive galaxies), since our classification is not in AGN 
luminosity cut, but more a cut in Lagn/^sf ratio. The recent result of 
Mullaney et al. (2012), that L X /L IR (which is a proxy of Lagn/^sf) 
is on average almost independent of stellar mass, implies that our 
selected AGN are not necessarily hosted by very massive galaxies, 
as indeed we find. 

Our results seem to confirm those of Santini et al. (2012), finding 
evidence of a higher average SF activity in AGN hosts with respect 
to inactive galaxies and a more pronounced level of SF enhance- 
ment in the hosts of luminous AGN (i.e. our AGN1 and AGN2). 


Most of the starburst galaxies (i.e. with SEDs fitted by local “star- 
burst” templates) are classified on-MS at any redshifts, though they 
are always in the region between the MS and the threshold. The 
starburst population seems to enhance its sSFR from z ~ 2 to z 
~ 1, occupying the region around the MS at higher z and shift- 
ing to higher SFRs (peaking at approximately two times the MS) 
at lower z. The bulk of the spiral population remains around the 
MS at all redshifts, although they almost disappear from our sam- 
ple at z ~ 2. From Fig. 15 we also note that while the SF-AGN 
population as a whole occupies both the loci of MS and off-MS 
sources, when divided into SF-AGN(SB) and SF-AGN(Spiral), it 
shows a very clear segregation (suggesting a different mode of SF 
for the two sub-classes). In fact, SF-AGN(SB) are mainly concen- 
trated in the off-MS region, while the bulk of the SF-AGN(Spiral) is 
on-MS. 

If we mix together all the ‘ingredients’ presented above, we can 
try to give a global interpretation to our results in terms of galaxy 
evolution. In agreement with other Herschel findings (i.e. Shao 
et al. 2010; Lapi et al. 201 1 ; Mullaney et al. 2012; Page et al. 2012; 
Rosario et al. 2012; Santini et al. 2012), we propose the following 
twofold evolutionary scenario for galaxies and AGN (sketched in a 
cartoon in Fig. 20). 
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Figure 20. A cartoon showing a possible evolutionary scenario involving AGN and star-forming galaxies and leading to the formation, on the one hand, of 
local elliptical galaxies, on the other hand, of local spiral galaxies. Top: a strong starburst with a growing BH inside [SF-AGN(SB)] transforms into an AGN 
(either AGN1 or AGN2, depending on orientation) when the BH mass and AGN luminosity have grown enough and the SF is likely to be quenched by feedback 
processes. Thereafter, the galaxy evolves passively towards a local elliptical. Bottom: the initial moderate starburst (starburst, lasting a typical time of ~few 
10 8 yr) transforms in a less intense starburst also fuelling a low-luminosity AGN [SF-AGN(Spiral)], which starts to reveal itself when the starburst activity 
is fading. Once the AGN is triggered, it heats the dust resulting in an increase in the mid-IR luminosity (flattening of the SED between 3 and 8 pm). The 
SF-AGN(Spiral) system could last about seven times longer than the pure starburst phase before becoming a steady spiral galaxy (spiral). 


(i) On the one hand, we observe that the sources with AGN- 
dominated SEDs (either AGN 1 or AGN2) and those with a starburst- 
like SED, but containing a non-dominant AGN [SF-AGN(SB)] have 
a peak in number and luminosity density at z ~ 2-2.5, dominating 
at higher redshifts and rapidly decreasing at lower z. The evolution 
of AGN1 and AGN2, both in luminosity and in density, is very 
similar, suggesting the same nature for type 1 and 2 AGN, and the 
unbiased power of observations in the far-IR band towards orien- 
tation/obscuration (affecting both optical and X-ray observations). 
While dominating the mid-IR part of the SED of these objects, 
the AGN is not able to explain the high observed far-IR emission, 
which is mostly powered by SF (as for the SF-AGN(SB) popula- 
tion, where the AGN never dominates the energetic output, probably 
due to dust-obscuration). The hosts of these AGN appear to form 
stars in a very efficient way, placing a large fraction of them above 
the known stellar mass-SFR MS (see Fig. 15). AGN1, AGN2 and 
SF-AGN(SB) are likely the progenitors of the elliptical galaxies ob- 
served nowadays in optical and near-IR surveys, forming through 
an intense burst of SF occurring during major mergers or in dense 
nuclear star-forming regions (Granato et al. 2001 ; Daddi et al. 2010; 
Wuyts et al. 201 1), and then followed by a phase of nuclear activity 
during which their Super Massive Black Holes (SMBHs) grow (i.e. 
Hopkins et al. 2008a, b; Lapi et al. 2011). It is generally agreed 
that the SMBHs and their host galaxies are tightly related, with 
major mergers being considered the likely process responsible for 


transporting large amount of gas towards the centre of the merging 
system, feeding the SMBH and triggering the SF activity. After the 
intense starburst phase, the AGN are believed to suppress the SF 
(i.e. Di Matteo, Springel & Hernquist 2005), so that the remnant 
quickly evolves to a red massive spheroid. This picture is strongly 
supported by the recent Herschel results from PEP (Rosario et al. 
2012), HerMES (Page et al. 2012; Harrison et al. 2012) and H- 
ATLAS (Lapi et al. 2011) surveys. The latter work, in particular, 
has shown that the bright-end of the IR LFs and counts at high 
redshift (>1.5) are consistent with the picture (e.g. Granato et al. 
2001) predicting the presence of a population of strongly obscured, 
star-forming galaxies with SED appreciably different from those of 
the local starbursts. X-ray and mid-IR spectroscopic observations 
(e.g. Alexander et al. 2005, Mullaney et al. 2012) of sub-mm se- 
lected sources have revealed in many of them the presence of a 
growing SMBH, powering an obscured AGN. In this framework, 
the most likely evolutionary path envisages first an SF-AGN(SB) 
phase (star-forming galaxy with a growing SMBH inside), then an 
AGN-dominated phase, and finally the formation of an elliptical 
galaxy in passive evolution. 

(ii) On the other hand, we observe that a significant fraction of 
our IR- selected sources is constituted by moderately star-forming 
galaxies characterized by an SED similar to that of spiral galaxies, 
but also suggesting the presence of a low-luminosity AGN [SF- 
AGN(Spiral)], best fitted by local Seyfert 1.8/2 templates. The bulk 
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of these objects and their principal contribution to pi R are at interme- 
diate redshifts (1 < z < 2), while they decrease between z ~ 1 and 
z = 0, as the spiral population rises. Most of the SF-AGN(Spiral)-, 
spiral- and starburst-SED galaxies occupy the region around the 
SFR-stellar mass MS at any redshifts, suggesting a steady mode of 
SF rather than a ‘starburst’ one for these three populations, whose 
evolution is likely connected. Indeed, Mullaney et al. (2012) have 
recently shown that at least up to z ~ 2, SMBHs have grown together 
with their host galaxies in star-forming galaxies, irrespective of host 
galaxy mass and triggering mechanism. Given the number densities 
(Fig. 18), evolutionary trends (Fig. 11), SFRs and masses (Figs 15 
and 19) of the SF-AGN(Spiral), starburst and spiral populations, we 
suggest that these three classes of objects might constitute different 
phases in the life of a galaxy undergoing secular evolution. The gas 
in moderate starburst galaxies, undergoing a burst of enhanced SF 
due either to gravitational interactions or disc instabilities (typical 
burst duration time of the order of a few 10 (i) * * * * * * 8 yr), might also fuel a 
low-luminosity AGN, which starts to reveal itself when the starburst 
activity is fading. Given the relatively high stellar masses found for 
the bulk of the SF-AGN(Spiral) [log(Af/M 0 ) ~ 10-11, see Fig. 15], 
and the almost constant M B h/M* ratio (of ~l-2 x 10~ 3 ) recently 
suggested by Mullaney et al. (2012) for all the 0.5 < z < 2.5 star- 
forming galaxies, they are likely to contain relatively massive BHs 
(M B h ~ 10 7 -10 8 Mq). They can therefore have either low values 
of their radiative efficiency or low values of their accretion mass 
rate m (m = M/M ^ dd < 0.01). AGN with low radiative efficiency 
or low accretion mass rates are generally called radiatively ineffi- 
cient accretion flows (Narayan et al. 1995), which include also the 
advection-dominated accretion flow (Narayan 1994). Low m AGN 
are more difficult to detect, since they often are less luminous than 
their hosts. Because of this, the nuclear emission is diluted by the 
host galaxy’s emission and many of them are likely to be classi- 
fied as ‘normal galaxies’ in most surveys if the AGN luminosity 
is less than that of the host (e.g. Hopkins et al. 2009). Given the 
relative number densities of the starburst and SF-AGN(Spiral) pop- 
ulations at 1 < z < 2, we hypothesize a typical duration of the 
SF-AGN(Spiral) phase about seven times longer than the starburst 
one (typical burst duration ~10 8 yr). Then, after a typical time of 
~7 x 10 8 yr, the AGN activity stops and these objects, whose num- 
ber density decreases at z < 1 , are likely to become steady spiral 
galaxies (rapidly increasing between z ~ 1 and z = 0) at lower 
redshifts. 


6 CONCLUSIONS 

We have used the 70- , 100-, 160-, 25 0-, 350- and 5 00- jam data 
from the cosmological guaranteed time Herschel surveys, PEP and 
HerMES, in the GOODS-S and -N, ECDFS and COSMOS, to char- 
acterize the evolution of the IR LF and luminosity density of PACS- 
selected sources across the redshift range 0 < z < 4. Evolution 
is well constrained by our data up to z ~ 3, and strong hints of 
evolution are derived at 3 < z < 4. In this work, we have: 

(i) completely characterized the multi wavelength SEDs of the 

PEP sources by performing a detailed SED-fitting analysis and a 

comparison with known template library of IR populations. Sources 

have been classified, based on their broad-band SEDs, in five main 

classes: spiral, starburst, SF-AGN, AGN1 and AGN2. 

(ii) computed the rest-frame LFs at 35, 60 and 90 pan up to z ~ 

4 from the 70-, 100- and 160- jam selected samples, respectively. 

(iii) integrated the SEDs over A, rest = 8-1000 jam and computed 
the total IR LF up to z ~ 4 and studied its evolution with redshift, 


finding strong luminosity evolution oc (1 + z ) 3 - 55±01 ° up to z ~ 
1.85, and a (1 + Z ) L62±0 - 51 between z ~ 1.85 andz ~ 4, combined 
with a negative density evolution oc (1 + z) -0 57 31 0-22 up to z ~ LI 
and oc (1 + z) -3 92 ± 0 34 at z > 1.1 and up to z ~ 4. 

(iv) derived the evolution of the comoving total IR luminosity 
density, which is found to increase as (1 + z) 30 ± 0-2 up to z ~ LI, 
then to remain nearly constant [decrease as (1 + z )-°- 3±01 ] from 
z ~ 1.1 to z ~ 2.8, and to decrease as (1 + z) -6 0 ± 0-9 up to z ~ 4. 

(v) found that the evolution derived for the global IR LF is indeed 
a combination of different evolutionary paths: the IR population 
does not evolve all together ‘as a whole’, as is often assumed in 
the literature, but is composed of different galaxy classes evolving 
differently: the spiral-SED galaxies dominate pi R only at low red- 
shifts (z < 0.5-0.6), then SF-AGN dominate up to z ~ 2.5, while 
AGN 1 and AGN2 start dominating the IR luminosity density only at 
z > 2.5. 

(vi) derived the relative contribution to pi R of MS and off-MS 
sources, which keep nearly constant between z ~ 0.8 and z ~ 2.2, 
with the MS population always dominating. The contribution to pi R 
of the off-MS sources shows no significant signs of increase with z 
(from ~ 19 per cent at z ~ 0.8-1.25 to >22 per cent at z ~ 1. 8-2.2). 

(vii) derived very different evolutionary behaviour in terms of 
different contributions to pi R , for galaxies with different masses, 
reflecting the downsizing scenario (pi R peaks at higher redshift with 
increasing mass). Intermediate-mass objects [log(M/MQ) = 10-1 1] 
always dominate the IR luminosity density, increasing with redshift 
up to z ~ 1 , then remaining nearly constant at higher redshifts (at 
least up to z ~ 2.8), while the contribution of most massive objects 
increases even more rapidly with z (at z ~ 2 it was approximately 
five times higher than today) and continues to grow up to z ~ 3. 

(viii) described a possible twofold evolutionary scenario for IR 
sources: on the one hand, AGN1 and AGN2, representing the same 
population, after an intense starburst phase (due to a major merging 
event, appearing as SF-AGN(SB) SED galaxies), suppress the SF 
(and shine as AGN) and evolve to red massive spheroids; on the 
other hand, the SF-AGN(Spiral) galaxies represent a phase in the 
life of a star-forming galaxy, following a moderate burst of SF 
(starburst, with SEDs like those of local starburst galaxies) and 
preceding the formation of a steady spiral galaxy as we observed in 
the local Universe. 
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